In this paper, the entire downhole fluid sucker rod pump system is replaced with a viscoelastic vibration model, namely a third-order differential equation with an inhomogeneous forcing term. Both Kelvin and Maxwell viscoelastic models can be implemented along with the dynamic behaviors of a mass point attached to the viscoelastic model. By employing the time-dependent polished rod force measured with a dynamometer as the input to the viscoelastic dynamic model, we have obtained the displacement responses, which match closely with the experimental measurements in actual operations, through an iterative process. The key discovery of this work is the feasibility of the so-called inverse optimization procedure which can be utilized to identify the equivalent scaling factor and viscoelastic system parameters. The proposed Newton-Raphson iterative method with some terms in the Jacobian matrix expressed with averaged rates of changes based on perturbations of up to two independent parameters provides a feasible tool for optimization issues related to complex engineering problems with mere information of input and output data from either experiments or comprehensive simulations. The same inverse optimization procedure is also implemented to model the entire fluid delivery system of a very viscous non-Newtonian polymer modeled as a first-order ordinary differential system similar to the transient entrance developing flow. The convergent parameter reproduces transient solutions which match very well with those from full-fledged computational fluid dynamics models with the required inlet volume flow rate and the outlet pressure conditions.