TOMLAB Solving Unconstrained and Constrained Optimization Problems
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. See 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,
Failed to parse (unknown function "\label"): {\displaystyle \label{rosen} \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 a = 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 6.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 (17) assuming a = 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 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 (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 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 (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 = 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 154 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 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 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 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 a 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 a 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 (17) 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,
Failed to parse (unknown function "\label"): {\displaystyle \label{ExpIII} \begin{array}{ll} \min\limits_{x} & f(x)= exp(x_1) (4 x_1^2+ 2 x_2^2+ 4 x_1 x_2+2 x_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 Section 7. 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 s 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 s = 0.2). The test for the last two cases are made for s = 0.9. If the user really wishes these methods to use s = 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
