HS071 Optimization Problem

For a solution to the HS071 constrained function minimization problem using a Python script instead of a notebook, see the HS071 Python script documentation.

Description

HS071 is Problem 71 from the collection of nonlinear programming test problems by Hoch and Schittkowski [1], who cite Bartholomew-Biggs [2] as the original source. The problem is to minimize the function

\[f(x)=x_{0} x_{3}\left(x_{0}+x_{1}+x_{2}\right)+x_{2}\]

subject to the nonlinear constraints

\[\begin{split}\begin{align} x_{0} x_{1} x_{2} x_{3} &\geq 25 \\ x_{0}^{2}+x_{1}^{2}+x_{2}^{2}+x_{3}^{2} & = 40 \end{align}\end{split}\]

and the variable bounds

\[1 \leq x_{i} \leq 5, \quad i=0, 1, 2, 3\]

The initial guess is given by

\[\boldsymbol{x} = [1,5,5,1]^T\]

(Note that we have used 0-based indexing in the problem statement, whereas Hoch and Schittkowski used 1-based indexing.)

YAPSS Solution

First, we instantiate the problem. It has no phases, 4 parameters, and 2 discrete constraints:

[1]:
from yapss import Problem

problem = Problem(name="HS071", nx=[], ns=4, nd=2)

We then define the objective and discrete constraint callback functions:

[2]:
def objective(arg):
    x = arg.parameter
    arg.objective = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2]


def discrete(arg):
    x = arg.parameter
    arg.discrete[:] = (
        x[0] * x[1] * x[2] * x[3],
        x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3],
    )


problem.functions.objective = objective
problem.functions.discrete = discrete

Set the bounds on the parameters and the constraint functions, per the problem statement:

[3]:
problem.bounds.parameter.lower = [1.0, 1.0, 1.0, 1.0]
problem.bounds.parameter.upper = [5.0, 5.0, 5.0, 5.0]
problem.bounds.discrete.lower = [25.0, 40.0]
problem.bounds.discrete.upper[1] = 40.0

We also provide an initial guess for the parameter values:

[4]:
problem.guess.parameter = [1.0, 5.0, 5.0, 1.0]

Specify YAPSS options, and also Ipopt options:

[5]:
problem.derivatives.order = "second"
problem.derivatives.method = "auto"
problem.ipopt_options.print_user_options = "no"
problem.ipopt_options.print_level = 3

Solve the problem:

[6]:
solution = problem.solve()

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

Total number of variables............................:        4
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        4
                     variables with only upper bounds:        0
Total number of equality constraints.................:        1
Total number of inequality constraints...............:        1
        inequality constraints with only lower bounds:        1
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0


Number of Iterations....: 8

                                   (scaled)                 (unscaled)
Objective...............:   1.7014017140224134e+01    1.7014017140224134e+01
Dual infeasibility......:   2.2928101314633036e-14    2.2928101314633036e-14
Constraint violation....:   2.1316282072803006e-14    2.1316282072803006e-14
Variable bound violation:   9.9907857542547163e-09    9.9907857542547163e-09
Complementarity.........:   1.0023967333275279e-11    1.0023967333275279e-11
Overall NLP error.......:   1.0023967333275279e-11    1.0023967333275279e-11


Number of objective function evaluations             = 9
Number of objective gradient evaluations             = 9
Number of equality constraint evaluations            = 9
Number of inequality constraint evaluations          = 9
Number of equality constraint Jacobian evaluations   = 9
Number of inequality constraint Jacobian evaluations = 9
Number of Lagrangian Hessian evaluations             = 8
Total seconds in IPOPT (w/o function evaluations)    =      0.013
Total seconds in NLP function evaluations            =      0.007

EXIT: Optimal Solution Found.

Finally, we format and print the solution:

[7]:
def print_variable(name, values):
    for i, _value in enumerate(values):
        print(f"{name}[{str(i)}] = {_value:1.6e}")


x = solution.parameter
print()

print("Solution of the primal variables, x")
print_variable("x", x)
print("\nSolution of the bound multipliers, z_L and z_U")
nlp_info = solution.nlp_info
print_variable("z_L", nlp_info.mult_x_L)
print_variable("z_U", nlp_info.mult_x_U)
print("\nSolution of the constraint multipliers, lambda")
print_variable("lambda", solution.discrete_multiplier)
print("\nObjective value")
print(f"f(x*) = {solution.objective:1.6e}")

Solution of the primal variables, x
x[0] = 1.000000e+00
x[1] = 4.743000e+00
x[2] = 3.821150e+00
x[3] = 1.379408e+00

Solution of the bound multipliers, z_L and z_U
z_L[0] = 1.087871e+00
z_L[1] = 2.671608e-12
z_L[2] = 3.544911e-12
z_L[3] = 2.635449e-11
z_U[0] = 2.499943e-12
z_U[1] = 3.896984e-11
z_U[2] = 8.482036e-12
z_U[3] = 2.762163e-12

Solution of the constraint multipliers, lambda
lambda[0] = -5.522937e-01
lambda[1] = 1.614686e-01

Objective value
f(x*) = 1.701402e+01

The solution above replicates the example given (in C++) in the Ipopt interface documentation, in much less code.

References