Solving Unconstrained and Constrained Optimization Problems

6 Solving Unconstrained and Constrained Optimization Problems

This section describes how to define and solve unconstrained and constrained optimization problems. Several examples are given on how to proceed, depending on if a quick solution is wanted, or more advanced runs are needed. TOMLAB is also compatible with MathWorks Optimization TB. See Appendix E for more information and test examples.

All demonstration examples that are using the Tomlab (TQ) format are collected in the directory examples . Running the menu program tomMenu , it is possible to run all demonstration examples. It is also possible to run each example separately. The examples relevant to this section are ucDemo and conDemo . All files that show how to use the Init File format are collected in the directory usersguide . The full path to these files are always given in the text. Throughout this section the test problem Rosenbrock's banana,

min

f(x)=α

x₂−x₁²

²+

1−x₁

s/t

−10	≤	x₁	≤	2
−10	≤	x₂	≤	2

(17)

is used to illustrate the solution of unconstrained problems. The standard value is α=100. In this formulation simple bounds are added on the variables, and also constraints in illustrative purpose. This problem is called RB BANANA in the following descriptions to avoid mixing it up with problems already defined in the problem definition files.

6.1 Defining the Problem in Matlab m-files

TOMLAB demands that the general nonlinear problem is defined in Matlab m-files, and not given as an input text string. A file defining the function to be optimized must always be supplied. For linear constraints the constraint coefficient matrix and the right hand side vector are given directly. Nonlinear constraints are defined in a separate file. First order derivatives and second order derivatives, if available, are stored in separate files, both function derivatives and constraint derivatives.

TOMLAB is compatible with MathWorks Optimization TB, which in various ways demands both functions, derivatives and constraints to be returned by the same function. TOMLAB handle all this by use of interface routines, hidden for the user. The user must then always use the MathWorks Optimization TB type of calls, not the TOMLAB function calls, and access to the GUI, menu and driver routines are not possible.

It is generally recommended to use the TOMLAB format instead, because having defined the files in this format, all MathWorks Optimization TB solvers are accessible through the TOMLAB multi-solver driver routines.

The rest of this section shows how to make the m-files for the cases of unconstrained and constrained optimization. These files does not depend on if the TQ or IF format are used to solve the problem, in both cases they are identical. The m-files for a constrained IF format example is shown. The files are defined in the directory usersguide and described in more detail in Appendix D. In Section 6.2 and onwards similar m-files are used to solve unconstrained optimization using the TQ format.

The most simple way to write the m-file to compute the function value is shown for the example in (17) assuming α=100:

File: tomlab/usersguide/rbbs_f.m

% crbb_f - function value for Constrained Rosenbrocks Banana
%
% function f = crbb_f(x)

function f = crbb_f(x)

f = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;

Running TOMLAB it is recommended to use a more general format for the m-files, adding one extra parameter, the Prob problem definition structure described in detail in Appendix A. TOMLAB will handle the simpler format also, but the more advanced features of TOMLAB are not possible to use.

If using this extra parameter, then any information needed in the low-level function computation routine may be sent as fields in this structure. For single parameter values, like the above α parameter in the example, the field Prob.user is recommended. See the Section 14.2 on how the to set User Supplied Problem Parameters in the Init File format in TOMLAB .

Using the above convention, then the new m-file for the example in (17) is defined as

File: tomlab/usersguide/rbb_f.m

% rbb_f - function value for Rosenbrocks Banana, Problem RB BANANA
%
% function f = rbb_f(x, Prob)

function f = rbb_f(x, Prob)

alpha = Prob.user.alpha;

f = alpha*(x(2)-x(1)^2)^2 + (1-x(1))^2;

The value in the field Prob.user is the α value. It is defined before calling the solver, either in a TOMLAB Init File, or directly before the call by explicit setting the Prob structure. In a similar way the gradient routine is defined as

File: tomlab/usersguide/rbb_g.m

% rbb_g - gradient vector for Rosenbrocks Banana, Problem RB BANANA
%
% function g = rbb_g(x, Prob)

function g = rbb_g(x, Prob)

alpha = Prob.user.alpha;

g = [-4*alpha*x(1)*(x(2)-x(1)^2)-2*(1-x(1)); 2*alpha*(x(2)-x(1)^2)];

If the gradient routine is not supplied, TOMLAB will use finite differences (or automatic differentiation) if the gradient vector is needed for the particular solver. In this case it is also easy to compute the Hessian matrix, which gives the following code

File: tomlab/usersguide/rbb_H.m

% rbb_H - Hessian matrix for Rosenbrocks Banana, Problem RB BANANA
%
% function H = crbb_H(x, Prob)

function H = crbb_H(x, Prob)

alpha = Prob.user.alpha;

H = [ 12*alpha*x(1)^2-4*alpha*x(2)+2 , -4*alpha*x(1);
                -4*alpha*x(1)        ,    2*alpha  ];

If the Hessian routine is not supplied, TOMLAB will use finite differences (or automatic differentiation) if the Hessian matrix is needed for the particular solver. Often a positive-definite approximation of the Hessian matrix is estimated during the optimization , and the second derivative routine is then not used.

If using the constraints defined for the example in (17) then a constraint routine needs to defined for the single nonlinear constraint, in this case

File: tomlab/usersguide/rbb_c.m

% rbb_c - nonlinear constraint vector for Rosenbrocks Banana, Problem RB BANANA
%
% function c = crbb_c(x, Prob)

function c = crbb_c(x, Prob)

cx = -x(1)^2 - x(2);

The constraint Jacobian matrix is also of interest and is defined as

File: tomlab/usersguide/rbb_dc.m

% rbb_dc - nonlinear constraint gradient matrix
%          for Rosenbrocks Banana, Problem RB BANANA
%
% function dc = crbb_dc(x, Prob)

function dc = crbb_dc(x, Prob)

% One row for each constraint, one column for each variable.

dc = [-2*x(1),-1];

If the constraint Jacobian routine is not supplied, TOMLAB will use finite differences (or automatic differentiation) to estimate the constraint Jacobian matrix if it is needed for the particular solver.

The solver nlpSolve is also using the second derivatives of the constraint vector. The result is returned as a weighted sum of the second derivative matrices with respect to each constraint vector element, the weights being the Lagrange multipliers supplied as input to the routine. For the example problem the routine is defined as

File: tomlab/usersguide/rbb_d2c.m

% rbb_d2c - The second part of the Hessian to the Lagrangian function for the
%     nonlinear constraints for Rosenbrocks Banana, Problem RB BANANA,i.e.
%
%           lam' * d2c(x)
%
% in
%
%   L(x,lam) =   f(x) - lam' * c(x)
% d2L(x,lam) = d2f(x) - lam' * d2c(x) = H(x) - lam' * d2c(x)
%
% function d2c=crbb_d2c(x, lam, Prob)

function d2c=crbb_d2c(x, lam, Prob)

% The only nonzero element in the second derivative matrix for the single
% constraint is the (1,1) element, which is a constant -2.

d2c = lam(1)*[-2 0; 0 0];

6.1.1 Communication between user routines

It is often the case that mathematical expressions that occur in the function computation also is part of the gradient and Hessian computation. If these operations are costly it is natural to avoid recomputing these and reuse them when computing the gradient and Hessian.

The function routine is always called before the gradient routine in TOMLAB , and the gradient routine is always called before the Hessian routine. The constraint routine is similarly called before the computation of the constraint gradient matrix. However, the TOM solvers call the function before the constraint routine, but the SOL solvers do the reverse.

Thus it is safe to use global variables to communicate information from the function routine to the gradient and Hessian, and similarly from the constraint routine to the constraint gradient routine. Any non-conflicting name could be used as global variable, see Table 30 in Appendix C to find out which names are in use. However, the recommendation is to always use a predefined global variable named US_A for this communication. TOMLAB is designed to handle recursive calls, and any use of new global variables may cause conflicts. The variable US_A (and also US_B ) is automatically saved in a stack, and any level of recursions may safely be used. The user is free to use US_A both as variable, and as a structure. If much information is to be communicated, defining US_A as a structure makes it possible to send any amount of information between the user routines.

In the examples directory the constrained optimization example in condemo is using the defined functions con1_f , con1_g and con1_H . They include an example of communicating one exponential expression between the routines.

The lsdemo example file in the examples directory communicates two exponential expressions between ls1_r and ls1_J with the use of US_A and US_B . In ls1_r the main part is

...
global US_A

t  = Prob.LS.t(:);

% Exponential computations takes time, and may be done once, and
% reused when computing the Jacobian
US_A = exp(-x(1)*t);
US_B = exp(-x(2)*t);

r = K*x(1)*(US_B - US_A) / (x(3)*(x(1)-x(2))) - Prob.LS.y;

In ls1_J then US_A is used

...
global US_A
% Pick up the globally saved exponential computations
e1 = US_A;
e2 = US_B;

% Compute the three columns in the Jacobian, one for each of variable
J = a * [ t.*e1+(e2-e1)*(1-1/b), -t.*e2+(e2-e1)/b, (e1-e2)/x(3)];

For more discussions on global variables and the use of recursive calls in TOMLAB , see Appendix C.

In the following sections it is described how to setup problems in TOMLAB and use the defined m-files. First comes the simplest way, to use the TOMLAB format.

6.2 Unconstrained Optimization Problems

The use of the TOMLAB format is best illustrated by examples

The following is the first example in the ucDemo demonstration file. It shows an example of making a call to probAssign to create a structure in the TOMLAB TQ format, and solve the problem with a call to ucSolve .


% ---------------------------------------------------------------------
function uc1Demo
% ---------------------------------------------------------------------

format compact
fprintf('=====================================================\n');
fprintf('Rosenbrocks banana with explicit f(x), g(x) and H(x)\n');
fprintf('=====================================================\n');

Name    = 'RB Banana';
x_0     = [-1.2 1]';   % Starting values for the optimization.
x_L     = [-10;-10];   % Lower bounds for x.
x_U     = [2;2];       % Upper bounds for x.
fLowBnd = 0;           % Lower bound on function.

% Generate the problem structure using the TOMLAB format (short call)
Prob    = probAssign('uc', x_L, x_U, Name, x_0, fLowBnd);

% Update the Prob structure with the names of files
Prob    = tomFiles(Prob,'uc1_f', 'uc1_g', 'uc1_H');

Result  = tomRun('ucSolve', Prob, 1);

Instead of using probAssign and tomFiles , it is possible to use conAssign with a limited number of input parameters. In its more general form conAssign is used to define constrained problems. It also takes as input the nonzero pattern of the Hessian matrix, stored in the matrix HessPattern . In this case all elements of the Hessian matrix are nonzero, and either HessPattern is set as empty or as a matrix with all ones. Also the parameter pSepFunc should be set. It defines if the objective function is partially separable, see Section 14.5. Setting this parameter empty (the default), then this property is not used. In the above example the call would be

...
HessPattern = ones(2,2);   % The pattern of nonzeros
pSepFunc    = [];          % No partial separability used

% conAssign is used to generate the TQ problem structure
Prob    = conAssign('uc1_f', 'uc1_g', 'uc1_H', HessPattern, ...
                     x_L, x_U, Name, x_0, pSepFunc, fLowBnd);
...

Also see the other examples in ucDemo on how to solve the problem, when gradients routines are not present, and numerical differentiation must be used. An example on how to solve a sequence of problems is also presented.

If the gradient is not possible to define, it is just to set the corresponding gradient function name empty, or reduce the number of parameters in the call to tomFiles , as the following example (uc2Demo ) shows:

...
% Only give the function. TOMLAB then estimates any derivatives automatically
Prob   = tomFiles(Prob,'uc1_f');
...

The example uc3Demo in file ucDemo show how to solve a sequence of problems in TOMLAB, in this case changing the steepness parameter α in (17). It is important to point out that it is only necessary to define the Prob structure once and then just change the varying parameters, in this case the α value. The version below is slightly modified, doing the call to conAssign making the parameter definitions directly. The α value is sent to the user routines using the field user in the Prob structure. Any field in the Prob structure could be used that is not conflicting with the predefined fields. In this example the a vector of Result structures are saved for later preprocessing.

% ---------------------------------------------------------------------
function uc3Demo - compact, slightly modified, version
% ---------------------------------------------------------------------

% conAssign is used to generate the TQ problem structure
% Prob = conAssign(f,g,H, HessPattern, x_L, x_U, Name, x_0, pSepFunc, fLowBnd);

Prob = conAssign('uc3_f',[],[],[],[-10;-10], [2;2], [-1.2;1], 'RB Banana',[],0)

% The different steepness parameters to be tried
Steep = [100 500 1000 10000];

for i = 1:4
    Prob.user.alpha      = Steep(i);
    Result(i)            = tomRun('ucSolve', Prob, 1);
end

6.3 Direct Call to an Optimization Routine

When wanting to solve a problem by a direct call to an optimization routine there are two possible ways of doing it. The difference is in the way the problem dependent parameters are defined. The most natural way is to use a Init File, like the predefined TOMLAB Init Files ◇_prob (e.g. uc_prob if the problem is of the type unconstrained) to define those parameters. The other way is to define those parameters by first calling the routines ProbAssign and tomFiles , or the routine conAssign . In this subsection, examples of two different approaches are given.

First, solve the problem RB BANANA in (17) as an unconstrained problem. In this case, define the problem in the files ucnew_prob , ucnew_f , ucnew_g and ucnew_H as described in Appendix D.3. Using the problem definition files in the working directory solve the problem and print the result by the following calls.

File: tomlab/usersguide/ucnewSolve1.m

probFile = 'ucnew_prob';       % Problem definition file.
P = 18;                        % Problem number for the added RB Banana.
Prob = probInit(probFile, P);  % Setup Prob structure.

Result = tomRun('ucSolve', Prob, 1);

Now, solve the same problem as in the example above but define the parameters x_0 , x_L and x_L by calling the routine ProbAssign . Note that in this case the file ucnew_prob is not used, only the files ucnew_f and ucnew_g . The file ucnew_H is not needed because a quasi-Newton BFGS algorithm is used. The call to the routine tomFiles defines the files that defines the problem.

File: tomlab/usersguide/ucnewSolve2.m

oType   = 'uc';        % Problem type.
x_0     = [-1.2;1];    % Starting values for the optimization.
x_L     = [-10;-10];   % Lower bounds for x.
x_U     = [2;2];       % Upper bounds for x.

Prob    = probAssign(oType, x_L, x_U, 'ucNew',x_0);% Setup Prob structure.
Prob    = tomFiles(Prob,'ucnew_f','ucnew_g');      % Problem definition files.

Prob.P  = 18;                                      % Problem number.
Prob.Solver.Alg=2;                                 % Use quasi-Newton BFGS
Prob.user.alpha = 100;                             % Set alpha parameter

Result = tomRun('ucSolve', Prob, 1);

Note that the calls to ProbAssign and tomFiles could be replaced with the following call to conAssign .

Prob = conAssign('ucnew_f','ucnew_g',[],[],[-10;-10], [2;2], 'ucNew', [-1.2;1]);

6.4 Constrained Optimization Problems

Study the following constrained exponential problem, Exponential problem III,

min

f(x)= exp(x₁) (4 x₁²+ 2 x₂²+ 4 x₁ x₂+2 x₂+1)

s/t

−10	≤	x₁	≤	10
−10	≤	x₂	≤	10
0	≤	x₁ + x₂	≤	0
1.5	≤	−x₁ x₂+x₁+x₂
−10	≤	x₁ x₂

. (18)

The first two constraints are simple bounds, the third is a linear equality constraint, because lower and upper bounds are the same. The last two constraints are nonlinear inequality constraints. To solve the problem, define the following statements, available as con1Demo in file conDemo .

Name    = 'Exponential problem III';
A       = [1 1];        % One linear constraint
b_L     = 0;            % Lower bound on linear constraint
b_U     = 0;            % b_L == b_U implies equality
c_L     = [1.5;-10]     % Two nonlinear inequality constraints
c_U     = [];           % Empty means Inf (default) for the two constraints
x_0     = [-5;5];       % Initial value for x
x_L     = [-10;-10];    % Lower bounds on x
x_U     = [10;10];      % Upper bounds on x
fLowBnd = 0;            % A lower bound on the optimal function value
x_min   = [-2;-2];      % Used for plotting, lower bounds
x_max   = [4;4];        % Used for plotting, upper bounds

x_opt=[-3.162278, 3.162278; -1.224745, 1.224745]; % Two stationary points
f_opt=[1.156627; 1.8951];

HessPattern = [];  % All elements in Hessian are nonzero.
ConsPattern = [];  % All elements in the constraint Jacobian are nonzero.
pSepFunc    = [];  % The function f is not defined as separable

% Generate the problem structure using the TOMLAB format
Prob = conAssign('con1_f', 'con1_g', 'con1_H', HessPattern, x_L, x_U, ...
          Name, x_0, pSepFunc, fLowBnd, A, b_L, b_U,  'con1_c', 'con1_dc',...
          [], ConsPattern, c_L, c_U, x_min, x_max, f_opt, x_opt);

Result  = tomRun('conSolve',Prob);
PrintResult(Result);

The following example, con2Demo in file conDemo , illustrates numerical estimates of the gradient and constrained Jacobian matrix. Only the statements different from the previous example is given. Note that the gradient routine is not given at all, but the constraint Jacobian routine is given. Setting Prob.ConsDiff greater than zero will overrule the use of the constraint Jacobian routine. The solver conSolve is in this case called directly.

% Generate the problem structure using the TOMLAB format
Prob = conAssign('con1_f', [], [], HessPattern,  x_L, x_U, Name, x_0, ...
          pSepFunc, fLowBnd, A, b_L, b_U, 'con1_c', 'con1_dc', [], ...
          ConsPattern, c_L, c_U, x_min, x_max, f_opt, x_opt);

Prob.NumDiff  = 1;   % Use standard numerical differences
Prob.ConsDiff = 5;   % Use the complex variable method to estimate derivatives

Prob.Solver.Alg = 0; % Use default algorithm in conSolve

Result  = tomRun('conSolve', Prob, 1);

The third example, con3Demo in file conDemo , shows how to solve the same problem for a number of different initial values on x. The initial values are stored in the matrix X0, and in each loop step Prob.x_0 is set to one of the columns in X0. In a similar way any of the values in the Prob structure may be changed in a loop step, if e.g. the loop is part of a control loop. The Prob structure only needs to be defined once. The different initial points reveal that this problem is nasty, and that several points fulfill the convergence criteria. Only the statements different from the previous example is given. A different solver is called dependent on which TOMLAB version is used.

X0              = [ -1 -5  1 0 -5 ;
                     1  7 -1 0  5];

% Generate the problem structure using the TOMLAB format
Prob = conAssign('con1_f', 'con1_g', 'con1_H', HessPattern, x_L, x_U, Name, ...
          X0(:,1), pSepFunc, fLowBnd, A, b_L, b_U, 'con1_c', 'con1_dc',...
          [], ConsPattern, c_L, c_U, x_min, x_max, f_opt, x_opt);

Prob.Solver.Alg = 0;
TomV            = tomlabVersion;

for i = 1:size(X0,2)
    Prob.x_0 = X0(:,i);

    if TomV(1:1) ~= 'M'
       % Users of v3.0 may instead call MINOS (or SNOPT, NPSOL in v3.0 /SOL)
       Result   = tomRun('minos',Prob, 2);
    else
       Result   = tomRun('conSolve',Prob, 2);
    end
end

The constrained optimization solvers all have proven global convergence to a local minimum. If the problem is not convex, then it is always difficult to assure that a global minimum has been reached. One way to make it more likely that the global minimum is found is to optimize very many times with different initial values. The fifth example, con5Demo in file conDemo illustrates this approach by solving the exponential problem 50 times with randomly generated initial points.

If the number of variables are not that many, say fifteen, another approach is to use a global optimization solver like glcSolve to crunch the problem and search for the global minimum. If letting it run long enough, it is very likely to find the global minimum, but maybe not with high precision. To run glcSolve the problem must be box-bounded, and the advise is to try to squeeze the size of the box down as much as possible. The sixth example, con6Demo in file conDemo , illustrates a call to glcSolve . It is very simple to do this call if the problem has been defined in the TOMLAB format. The statements needed are the following

Prob.optParam.MaxFunc = 5000; % Define maximal number of function evaluations
Result                = tomRun('glcSolve',Prob,2);

A more clever approach, using warm starts and successive checks on the best function value obtained, is discussed in Section 7. It is also better to use glcAssign and not conAssign if the intension is to use global optimization.

6.5 Efficient use of the TOMLAB solvers

To follow the iterations in the TOMLAB Base Module solvers, it is useful to set the IterPrint parameter as true. This gives one line of information for each iteration. This parameter is part of the optParam subfield:

Prob.optParam.IterPrint = 1;

Note that ucSolve implements a whole set of methods for unconstrained optimization. If the user explicitly wants Newtons method to be used, utilizing second order information, then set

Prob.Solver.Alg=1;            % Use Newtons Method

But ucSolve will switch to the default BFGS method if no routine has been given to return the Hessian matrix. If the user still wants to run Newtons method, then the Hessian routine must be defined and return an empty Hessian. That triggers a numerical estimation of the Hessian. Do help ucSolve to see the different algorithmic options and other comments on how to run the solver.

Both ucSolve and conSolve use line search based methods. The parameter σ influences the accuracy of the line search each step. The default value is

Prob.LineParam.sigma = 0.9;   % Inaccurate line search

However, using the conjugate gradient methods in ucSolve , they benefit from a more accurate line search

Prob.LineParam.sigma = 0.01;   % Default accurate line search for C-G methods

as do quasi-Newton DFP methods (default σ = 0.2). The test for the last two cases are made for σ = 0.9. If the user really wishes these methods to use σ = 0.9, the value must be set slightly different to fool the test:

Prob.LineParam.sigma = 0.9001; % Avoid the default value for C-G methods

The choice of line search interpolation method is also important, a cubic or quadratic interpolation. The default is to use cubic interpolation.

Prob.LineParam.LineAlg = 1; % 0 = quadratic, 1 = cubic

« Previous « Start » Next »