Optimization with Linear Equality Constraints#

This GAUSS optimization example demonstrates the use of COMT to perform optimization with linear equality constraints.

Key example features#

The user-defined objective function fct for computing the Woods Function.
The first case uses comtControl linear equality members, c0.A and c0.B, to specify equality constraints.
The second case uses a user-defined function eqp in combination with the c0.eqProc member of the comtControl structure to specify the equality constraints.
The c0.bounds member of the comtControl stucture is used to set all parameter values between -10 and 10.

There are four linear equality constraints that are implemented in this example:

\[`x_1 + 3x_2 = 0`\]

\[`x_3 + x_4 - 2x_5 = 0`\]

\[`x_2 - x_5 = 0`\]

where \(x_1, x_2, \ldots, x_5\) are the parameters to be estimated.

Case One: Specifying equality constraints with control structure matrices#

In this case, two members of the comtControl are used to specify equality constraints. These set up the constraints based on the relationship:

\[`A(\theta) = B`\]

The c0.A matrix is used to specify the LHS of the equality constraints:

Each row represents one constraint. In other words, the matrix should have the same number of rows as there are constraints.
Each column reprents one parameter. In other words, the matrix should have the same number of columns as there are parameters.

The c0.B vector is used to specify the RHS of the equality constraints:

Each row represents one constraint. The matrix should have the same number of rows as there are constraints.
Always has one column.

new;
cls;
library comt;

// Nonlinear objective function
proc fct(x, ind);
    struct modelResults mm;

    if ind[1];
        // Calculate objective function value
        mm.function = (x[1] - x[2])^2 + (x[2] + x[3]- 2)^2 + (x[4] - 1)^2 + (x[5] - 1)^2;
    endif;

    // Return modelResults struct
    retp(mm);
endp;

// Declare 'c0' to be a comtControl struct
// and fill with default settings
struct comtControl c0;
c0 = comtControlCreate();

/*
** Specify equality constraints
*/
// The c0.A matrix represents represents the LHS
// of the equality constraints.
// Each row represents one constraint.
// Each column represents a parameter.
// The matrix should have:
//   - the same number of rows as constraints
//   - the same number of columns as parameters
c0.A = {1 3 0 0 0,
        0 0 1 1 -2,
        0 1 0 0 -1};

c0.B = 0;

// Turn off trust radius
c0.trustRadius = error(0);

// Bound all parameters to be between -10 and +10
c0.bounds = { -10 10 };

// Use Newton descent method
c0.algorithm = 3;

// Starting parameter values
x0 = { 2, 2, 2, 2, 2 };

// Declare 'out' to be a comtResults
// struct to hold optimization results
struct comtResults out;

// Minimize objective function
out = comt(&fct, x0, c0);

// Print optimization results
call comtPrt(out);

There code prints results to the Command Window.

Results#

The first portion of the results provide details about convergence and performance.

Return code    = 0
Function value = 4.09302
Convergence    : normal convergence

These results indicate that the optimization converged normally, with a return code of 0. Any return Code other than 0 would indicate some issue with the convergence. The exact meaning of the return code can be found in the comt() documentation.

The next section of the results reports the parameter estimates and the associated gradients.

Parameters  Estimates    Gradient
---------------------------------------------------------------------
x[1,1]        -0.7675     -2.0465
x[2,1]         0.2558     -0.1861
x[3,1]         0.6279     -2.2326
x[4,1]        -0.1162     -2.2325
x[5,1]         0.2558     -1.4884

The final section of the results reports the number of iterations and computation time.

Number of iterations    25
Minutes to convergence  0.00028

Case Two: Specifying equality constraints with user-defined equality function#

In this case, the same optimization is performed, using a procedure to define the equality constraints. Whenan equation is used to specify constraints, a pointer to the function must be provided using the c0.eqProc member of the comtControl structure.

new;
cls;
library comt;

// Nonlinear objective function
proc fct(x, ind);
    struct modelResults mm;

    if ind[1];
        // Calculate objective function value
        mm.function = (x[1] - x[2])^2 + (x[2] + x[3] - 2)^2 + (x[4] - 1)^2 + (x[5] - 1)^2;
    endif;

    // Return modelResults struct
    retp(mm);
endp;

// Procedure to compute equality constraints
// this must specify the constraint such that
// eqp(x) = 0
proc eqp(x);
    local result;

    // This will be returned and
    // it should be a vector of zeros
    // with the same number of rows as constraints
    result = zeros(3, 1);

    // Constraint one
    result[1] = x[1] + 3*x[2];

    // Constraint two
    result[2] = x[3] + x[4] - 2 * x[5];

    // Constraint three
    result[3] = x[2] - x[5];

    retp(result);
endp;

// Declare 'c0' to be a comtControl struct
// and fill with default settings
struct comtControl c0;
c0 = comtControlCreate();

// Turn off trust radius
c0.trustRadius = error(0);

// Assign pointer for equality procedure
c0.eqProc = &eqp;

// Bound all parameters to be between -10 and +10
c0.bounds = { -10 10 };

// Use Newton descent method
c0.algorithm = 3;

// Starting parameter values
x0 = { 2, 2, 2, 2, 2 };

// Declare 'out' to be a comtResults
// struct to hold optimization results
struct comtResults out;

// Minimize objective function
out = comt(&fct,x0,c0);

// Print optimization results
call comtPrt(out);

Results#

We can confirm that these results match those from the first case.

Return code    = 0
Function value = 4.09302
Convergence    : normal convergence

Parameters  Estimates    Gradient
---------------------------------------------------------------------
x[1,1]        -0.7674     -2.0465
x[2,1]         0.2558     -0.1861
x[3,1]         0.6279     -2.2326
x[4,1]        -0.1163     -2.2325
x[5,1]         0.2558     -1.4884


Number of iterations    63
Minutes to convergence  0.00050