TOMLAB Solving Unconstrained and Constrained Optimization Problems

This page is part of the TOMLAB Manual. See TOMLAB Manual.

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. See TOMLAB Appendix E for more information on your external solvers.

All demonstration examples that are using the TOMLAB format are collected in the directory examples. It is also possible to run each example separately. The examples relevant to this section are ucDemo and conDemo. The full path to these files are always given in the text. Throughout this section the test problem Rosenbrock's banana,

Equation: Rosenbrock's banana

${\begin{array}{ll}\min \limits _{x}&f(x)=\alpha \left(x_{2}-x_{1}^{2}\right)^{2}+\left(1-x_{1}\right)^{2}\\s/t&{\begin{array}{lcccl}-10&\leq &x_{1}&\leq &2\\-10&\leq &x_{2}&\leq &2\\\end{array}}\end{array}}$

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.

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.

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. In Section 1.2 and onwards similar m-files are used to solve unconstrained optimization using the TOMLAB format.

The most simple way to write the m-file to compute the function value is shown for the example in #Equation: Rosenbrock's banana assuming α = 100:

File: tomlab/usersguide/rbbs_f.m

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

function f = rbbs_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 TOMLAB 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 a parameter in the example, the field Prob.user is recommended.

Using the above convention, then the new m-file for the example in #Equation: Rosenbrock's banana 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 a value. It is defined before calling the solver 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  = rbb_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 #Equation: Rosenbrock's banana 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 = rbb_c(x,  Prob) 

function cx  = rbb_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 = rbb_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=rbb_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];

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 in Appendix D 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 TOMLAB Appendix D.

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.

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 conAssign to create a structure in the TOMLAB 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];	
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 = conAssign('uc1_f', 'uc1_g',  'uc1_H', [], x_L,  x_U, Name, ... 
                 x_0,  [], fLowBnd);

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

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 TOMLAB Special Notes and Features. 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  TOMLAB  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 one just sets the corresponding gradient function name empty.

The example uc3Demo in file ucDemo show how to solve a sequence of problems in TOMLAB, in this case changing the steepness parameter α. 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 a 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  - Sequence  of  Rosenbrocks banana functions
%   ---------------------------------------------------------------------

%   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

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 o prob (e.g. uc prob if the problem is of the type unconstrained) to define those parameters. The other way is to call the routine conAssign. In this subsection, examples of two different approaches are given.

First, solve the problem RB BANANA in #Equation: Rosenbrock's banana as an unconstrained problem. In this case, define the problem in the files ucnew prob, ucnew f, ucnew g and ucnew H. 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 conAssign. 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.

File: tomlab/usersguide/ucnewSolve2.m

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 = conAssign('ucnew_f','ucnew_g', [], [], x_L,  x_U,...
                 'ucNew', x_0);

Prob.P 	= 18; 	                   %   Problem number. 
Prob.Solver.Alg=1;	           %   Use quasi-Newton  BFGS 
Prob.user.uP = 100; 	           %   Set  alpha  parameter 
Result = tomRun('ucSolve',Prob,1);

Constrained Optimization Problems

Study the following constrained exponential problem, Exponential problem III,

${\begin{array}{ll}\min \limits _{x}&f(x)=exp(x_{1})(4x_{1}^{2}+2x_{2}^{2}+4x_{1}x_{2}+2x_{2}+1)\\s/t&{\begin{array}{lcccl}-10&\leq &x_{1}&\leq &10\\-10&\leq &x_{2}&\leq &10\\0&\leq &x_{1}+x_{2}&\leq &0\\1.5&\leq &-x_{1}x_{2}+x_{1}+x_{2}&&\\-10&\leq &x_{1}x_{2}&&\\\end{array}}\end{array}}.$

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 X_0, and in each loop step Prob.x_0 is set to one of the columns in X_0. 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 TOMLAB Solving Global Optimization Problems. It is also better to use glcAssign and not conAssign if the intension is to use global optimization.

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