Fits data to a nonlinear model (possibly with linear constraints) using the successive quadratic programming algorithm (applied to the sum of squared errors, sse = Σ(yi − f(xi; θ))2) and either a finite difference gradient or a user-supplied gradient.

The type double function is imsls_d_nonlinear_optimization.

Function imsls_f_nonlinear_optimization is based on M.J.D. Powell’s TOLMIN, which solves linearly constrained optimization problems, i.e., problems of the form min f(θ), θ ∈ ℜ, subject to

given the vectors b1,  b2, θI, and θu and the matrices A1 and A2.

The algorithm starts by checking the equality constraints for inconsistency and redundancy. If the equality constraints are consistent, the method will revise θ0, the initial guess provided by the user, to satisfy

Next, θ0 is adjusted to satisfy the simple bounds and inequality constraints. This is done by solving a sequence of quadratic programming subproblems to minimize the sum of the constraint or bound violations.

Now, for each iteration with a feasible θk, let Jk be the set of indices of inequality constraints that have small residuals. Here, the simple bounds are treated as inequality constraints. Let Ik be the set of indices of active constraints. The following quadratic programming problem

subject to

is solved to get (dk, λk) where aj is a row vector representing either a constraint in A1 or A2 or a bound constraint on θ. In the latter case, the aj = ei for the bound constraint θi ≤ (θu)i and aj = −ei for the constraint θi ≤ (θl)i. Here, ei is a vector with a 1 as the i‑th component, and zeroes elsewhere. λk are the Lagrange multipliers, and Bk is a positive definite approximation to the second derivative ∇2 f(θk).

After the search direction dk is obtained, a line search is performed to locate a better point. The new point θk+1 = θk + αkdk has to satisfy the conditions

and

The main idea in forming the set Jk is that, if any of the inequality constraints restricts the step-length αk, then its index is not in Jk. Therefore, small steps are likely to be avoided.

Finally, the second derivative approximation, Bk, is updated by the BFGS formula, if the condition

holds. Let θk ← θk+1, and start another iteration.

The iteration repeats until the stopping criterion

is satisfied; here, τ is a user-supplied tolerance. For more details, see Powell (1988, 1989).

Since a finite-difference method is used to estimate the gradient for some single precision calculations, an inaccurate estimate of the gradient may cause the algorithm to terminate at a noncritical point. In such cases, high precision arithmetic is recommended. Also, whenever the exact gradient can be easily provided, the gradient should be passed to imsls_f_nonlinear_optimization using the optional argument IMSLS_JACOBIAN.

In this example, a data set is fitted to the nonlinear model function

Draper and Smith (1981, p. 475) state a problem due to Smith and Dubey. [H. Smith and S. D. Dubey (1964), "Some reliability problems in the chemical industry", Industrial Quality Control, 21 (2), 1964, pp. 64−70] A certain product must have 50% available chlorine at the time of manufacture. When it reaches the customer 8 weeks later, the level of available chlorine has dropped to 49%. It was known that the level should stabilize at about 30%. To predict how long the chemical would last at the customer site, samples were analyzed at different times. It was postulated that the following nonlinear model should fit the data.

Since the chlorine level will stabilize at about 30%, the initial guess for theta1 is 0.30. Using the last data point (x = 42, y = 0.39) and θ0 = 0.30 and the above nonlinear equation, an estimate for θ1of 0.02 is obtained.

The constraints that θ0 ≥ = 0 and θ1 ≥ = 0 are also imposed. These are equivalent to requiring that the level of available chlorine always be positive and never increase with time.

The Jacobian of the nonlinear model equation is also used.