Getting Started

A basic example will help us to get started with WORHP. First of all, we need an optimisation problem to be solved.

Problem Formulation

Let the optimisation variables be $x \in \mathbb{R}^{4}$ . The objective function of interest is $f(x) := x_1^2 + 2 x_2^2 - x_3$ with respect to the nonlinear constraint $g_1(x) := x_1^2 + x_3^2 + x_1x_3$ and the linear constraints $g_2(x) := x_3 - x_4$ and $g_3(x) := x_2 + x_4$ . Furthermore, the optimisation variables must satisfy the box constraints $-0.5 \leq x_1,$ $-2 \leq x_2,$ $0 \leq x_3 \leq 2$ and $-2 \leq x_4 \leq 2$ .

WORHP Specific Formulation

$\begin{equation*} \begin{array}{rrcllr} & &\min\limits_{x\in \mathbb{R}^{4}}& x_1^2 + 2 x_2^2 - x_3&& \\ \\ \mbox{subject to}&& & x_1^2 + x_3^2 + x_1x_3 &=& 1 \\ &-\infty&\leq& x_3 - x_4 &\leq& -1 \\ &2.5 &\leq& x_2 + x_4 & \leq & 5 \\ &-0.5& \leq& x_1 & \leq & \infty \\ &-2& \leq& x_2 &\leq & \infty \\ &0 &\leq& x_3 & \leq & 2 \\ &-2& \leq& x_4 & \leq & 2 \end{array} \end{equation*}$

The WORHP specific formulation requires to mark constraints or variables without upper bound with infinity $(\infty)$ and without lower bounds with minus infinity $(-\infty)$ respectively.

Definition of Derivatives

WORHP uses derivative based optimisation to minimise the given problem. Thus, in order to improve the performance of the solver analytic derivatives should be given. In this section we explain the resulting derivatives of our example. The implementations in the different programming languages can be found in the next section. Remember to take care of the indexing of your programming language of choice (i.e. FORTRAN 1-based, C/C++ 0-based). We will use FORTRAN 1-based vectors within this section. Please have a look at the section "Sparse Matrices" within our manual for further information on sorting requirements for the matrix entries within WORHP. Additionally, WORHP uses a scaling factor for the objective function $\sigma$ (wsp%ScaleObj), which must be used when computing the objective function value, the value of the gradient of the objective and within the computation of the Hessian of the Lagrangian.

The gradient of the objective function $(\nabla f(x))$ for our example is given by $\begin{equation*} \nabla f(x) = \begin{pmatrix} 2x_1 \\ 4x_2 \\ -1 \\ 0 \end{pmatrix}. \end{equation*}$ The implementation uses sparse matrices for the derivatives. Therefore, this gradient leads to the row vector DFrow = [1, 2, 3] and the value vector DFval = [wsp%ScaleObj * 2.0 * x(1), wsp%ScaleObj * 4.0 * x(2), wsp%ScaleObj * (-1.0)]. If you do not provide an implementation of the objective gradient, you must set the parameter UserDF = False.

The jacobian of the constraint function $(\nabla g(x))$ is given by $\begin{equation*} \nabla g(x) = \begin{pmatrix} 2x_1+x_3 & 0 & 2x_3+x_1 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 1 & 0 & 1 \end{pmatrix}. \end{equation*}$ The row vector is DGrow = [1, 3, 1, 2, 2, 3], the column vector is DGcol = [1, 2, 3, 3, 4, 4] and the value vector is DGval = [2.0*x(1) + x(3), 1.0, 2.0*x(3)+x(1), 1.0, -1.0, 1.0]. If you do not provide an implementation of the jacobian of the constraints, you must set the parameter UserDG = False.

The implemented algorithms require the use of the Lagrangian function $L(x, \mu) := \sigma f(x) + \mu^\top g(x)$ , with the scaling factor $\sigma$ (wsp%ScaleObj) for the objective part within the Lagrangian. One advantage of the used Sequential Quadratic Programming and Interior-Point approach is, that theoretically local quadratic convergence is achievable. To reach this order of convergence second order derivatives of the Lagrangian must be used. Due to symmetry only the lower triangular part of the matrix must be stored. Additionally, WORHP requires a special sorting of the sparse entries. All entries on the diagonal must be given, even structural zeros. Furthermore, the non-diagonal entries must be given first. To get a deeper insight, please have a look at the WORHP manual. $\begin{equation*} \nabla^2_{xx} L(x, \mu) = \begin{pmatrix} 2\sigma + 2\mu_1 & 0 & \mu_1 & 0 \\ 0 & 4 \sigma& 0 & 0 \\ \mu_1 & 0 & 2\mu_1 & 0 \\ 0 & 0 &0 & 0 \end{pmatrix}. \end{equation*}$ The row vector is HMrow = [3, 1, 2, 3, 4], the column vector is HMcol = [1, 1, 2, 3, 4] and the value vector is HMval = [wsp%ScaleObj * Mu(1), 2.0 + 2.0*Mu(1), wsp%ScaleObj * 4.0, 0.0, 0.0]. If you do not provide an implementation of the Lagrangian of the Hessian, you must set the parameter UserHM = False.

Source Code

Please check out the source code of the language of your choice to get further insights:

C C++ FORTRAN Python

Output

After running the example, the output looks like this:

 ReadParamsNoInit: Used parameter file worhp.xml
 * Read 268/268 parameters.

 -------------------------------------------------------
     This is WORHP 1.13-0, the European NLP-solver.
    Use of WORHP is subject to terms and conditions.
     Visit http://www.worhp.de for more information.
 -------------------------------------------------------

 Total number of variables ........................        4
                                    fixed variables        0
                    variables with lower bound only        2
               variables with lower and upper bound        2
                    variables with upper bound only        0
 Total number of box constraints ..................        6
 Total number of other constraints ................        4
                               equality constraints        1
       inequality constraints with lower bound only        0
  inequality constraints with lower and upper bound        1
       inequality constraints with upper bound only        1


 Gradient (user)  3/4  =  75.000%
 Jacobian (user)  6/12 =  50.000%
 Hessian  (user)  5/10 =  50.000%

 Algorithm    Sequential Quadratic Programming
 NLP MaxIter  10000           Line Search Method     Filter
 QP  MaxIter  500
 LA solver    MA97 (tol 1.0E-09, ref 10, ord METIS/AMD, scl none)

 Tolerances:
 Optimality  (sKKT)   1.00E-06 (1.0E-03)   IP ComTol    2.00E-07
 Feasibility          1.00E-06 (1.0E-03)   IP ResTol    4.00E-08
 Complementarity      1.00E-06

 Timeout              1800.000  seconds

      ITER           OBJ             CON           OPTI/COMPL     FLAGS   ALPHA       |DX|     REL   PEN   REG     TIME
 [    0|   6]  +1.10000000E+01  6.00000000E+00  6.62989673E+00 sc  Uin  0.000E+00   2.77E+00    -     -     -    9.20E-04
 [    1|   3]  +1.38081440E-01  1.43688144E+00  8.28496615E-01 c   Uin  1.000E+00   2.77E+00    -     -   -2.90  1.25E-03
 [    2|   9]  -4.39931123E-01  3.05158405E-01  1.32548214E-01 c   Uin  1.000E+00   5.54E-01    -     -   -3.12  1.81E-03
 [    3|   1]  -4.98391574E-01  4.17134806E-02  2.19258787E-02 so  Uin  1.000E+00   2.05E-01    -     -   -3.45  2.01E-03
 [    4|   8]  -4.99997835E-01  1.47366536E-03  4.71936920E-04 so  Uia  1.000E+00   3.86E-02    -     -   -3.89  2.51E-03
 [    5|   2]  -5.00000000E-01  2.15886338E-06  6.49481731E-08 so  Uao  1.000E+00   1.47E-03    -     -   -4.32  2.76E-03
 [    6|   1]  -5.00000000E-01  4.96989561E-10  2.53618587E-15 so  Ufo  1.000E+00   2.16E-06   0.00   -   -4.59  2.97E-03


 Final values after iteration 6:
 Final objective value ............. -5.0000000000E-01
 Final constraint violation ........  4.9698956062E-10
 Final complementarity .............  0.0000000000E+00 (0.0000000000E+00)
 Final KKT conditions ..............  2.5361858705E-15 (1.1745104643E-09)
 Successful termination: Optimal Solution Found.

A detailed description of the output can be found in the Users' Guide.