TOMLAB Defining Problems in TOMLAB

From TomWiki
Jump to navigationJump to search

Notice.png

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

TOMLAB is based on the principle of creating a problem structure that defines the problem and includes all relevant information needed for the solution of the user problem. One unified format is defined, the TOMLAB format. The TOMLAB format gives the user a fast way to setup a problem structure and solve the problem from the Matlab command line using any suitable TOMLAB solver.

TOMLAB also includes a modeling engine (or advanced Matlab class), TomSym, see TOMLAB Defining Problems in TOMLAB#TomSym. The package uses Matlab objects and operator overloading to capture Matlab procedures, and generates source code for derivatives of any order.

In this section follows a more detailed description of the TOMLAB format.

The TOMLAB Format

The TOMLAB format is a quick way to setup a problem and easily solve it using any of the TOMLAB solvers. The principle is to put all information in a Matlab structure, which then is passed to the solver, which extracts the relevant information. The structure is passed to the user function routines for nonlinear problems, making it a convenient way to pass other types of information.

The solution process for the TOMLAB format has three steps:

  1. Define the problem structure, often called Prob.
  2. Call the solver or the universal driver routine tomRun.
  3. Postprocessing, e.g. print the result of the optimization.

Step 1 could be done in several ways in TOMLAB. Recommended is to call one of the following routines dependent on the type of optimization problem, see #Table: Routines to create a problem structure in the TOMLAB format..

Step 2, the solver call, is either a direct call, e.g. conSolve: 

Prob = ProbCheck(Prob,  'conSolve'); 
Result  = conSolve(Prob);

or a call to the multi-solver driver routine tomRun, e.g. for constrained optimization: 

Result  = tomRun('conSolve', Prob);

Note that tomRun handles several input formats. Step 3 could be a call to PrintResult.m: 

PrintResult(Result);

The 3rd step could be included in Step 2 by increasing the print level to 1, 2 or 3 in the call to the driver routine 

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

See the different demo files that gives examples of how to apply the TOMLAB format: conDemo.m, ucDemo.m, qpDemo.m, lsDemo.m, lpDemo.m, mipDemo.m, glbDemo.m and glcDemo.m.

Table: Routines to create a problem structure in the TOMLAB format.

Matlab call probTypes Type of optimization problem
Prob = bmiAssign( ... ) 14 Semi-definite programming with bilinear matrix inequalities.
Prob = clsAssign( ... ) 4,5,6 Unconstrained and constrained nonlinear least squares.
Prob = conAssign( ... ) 1,3 Unconstrained and constrained nonlinear optimization.
Prob = expAssign( ... ) 17 Exponential fitting problems.
Prob = glcAssign( ... ) 9,10,15 Box-bounded or mixed-integer constrained global programming.
Prob = lcpAssign( ... ) 22 Linear mixed-complimentary problems.
Prob = llsAssign( ... ) 5 Linear least-square problems.
Prob = lpAssign( ... ) 8 Linear programming.
Prob = lpconAssign( ... ) 3 Linear programming with nonlinear constraints.
Prob = mcpAssign( ... ) 23 Nonlinear mixed-complimentary problems.
Prob = minlpAssign( ... ) 12 Mixed-Integer nonlinear programming.
Prob = mipAssign( ... ) 7 Mixed-Integer programming.
Prob = miqpAssign( ... ) 11 Mixed-Integer quadratic programming.
Prob = miqqAssign( ... ) 18 Mixed-Integer quadratic programming with quadratic constraints.
Prob = qcpAssign( ... ) 23 Quadratic mixed-complimentary problems.
Prob = qpblockAssign( ... ) 2 Quadratic programming (factorized).
Prob = qpAssign( ... ) 2 Quadratic programming.
Prob = qpconAssign( ... ) 3 Quadratic programming with nonlinear constraints.
Prob = sdpAssign( ... ) 13 Semi-definite programming with linear matrix inequalities.
Prob = amplAssign( ... ) 1-3,7,8,11,12 For AMPL problems defined as nl -files.
Prob = simAssign( ... ) 1,3-6,9-10 General routine, functions and constraints calculated at the same time

Modifying existing problems

It is possible to modify an existing Prob structure by editing elements directly, however this is not recommended since there are dependencies for memory allocation and problem sizes that the user may not be aware of.

There are a set of routines developed specifically for modifying linear constraints (do not modify directly, Prob.mLin need to be set to a proper value if so). All the static information can be set with the following routines.

add_A

Purpose

Adds linear constraints to an existing problem.

Calling Syntax 

Prob = add_A(Prob, A, b_L, b_U)

Description of Inputs

Prob Existing TOMLAB problem.
A The additional linear constraints.
b_L The lower bounds for the new linear constraints.
b_U The upper bounds for the new linear constraints.

Description of Outputs

Prob Modified TOMLAB problem.

keep_A

Purpose

Keeps the linear constraints specified by idx.

Calling Syntax 

Prob = keep_A(Prob, idx)

Description of Inputs

Prob Existing TOMLAB problem.
idx The row indices to keep in the linear constraints.

Description of Outputs

Prob Modified TOMLAB problem.

remove A

Purpose

Removes the linear constraints specified by idx.

Calling Syntax 

Prob = remove_A(Prob, idx)

Description of Inputs

Prob Existing TOMLAB problem.
idx The row indices to remove in the linear constraints.

Description of Outputs

Prob Modified TOMLAB problem.

replace A

Purpose

Replaces the linear constraints.

Calling Syntax 

Prob = replace_A(Prob, A, b_L, b_U)

Description of Inputs

Prob Existing TOMLAB problem.
A New linear constraints.
b_L Lower bounds for linear constraints.
b_U Upper bounds for linear constraints.

Description of Outputs

Prob Modified TOMLAB problem.

modify b_L

Purpose

Modify lower bounds for linear constraints. If idx is not given b_L will be replaced.

Calling Syntax 

Prob = modify_b_L(Prob, b_L, idx)

Description of Inputs

Prob Existing TOMLAB problem.
b_L New lower bounds for the linear constraints.
idx Indices for the modified constraint bounds (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify b_U

Purpose

Modify upper bounds for linear constraints. If idx is not given b_U will be replaced.

Calling Syntax 

Prob = modify_b_U(Prob, b_U, idx)

Description of Inputs

Prob Existing TOMLAB problem.
b_U New upper bounds for the linear constraints.
idx Indices for the modified constraint bounds (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify c

Purpose

Modify linear objective (LP/QP only).

Calling Syntax 

Prob = modify_c(Prob, c, idx)

Description of Inputs

Prob Existing TOMLAB problem.
c New linear coefficients.
idx Indices for the modified linear coefficients (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify c_L

Purpose

Modify lower bounds for nonlinear constraints. If idx is not given c L will be replaced.

Calling Syntax 

Prob = modify_c_L(Prob, c_L, idx)

Description of Inputs

Prob Existing TOMLAB problem.
c_L New lower bounds for the nonlinear constraints.
idx Indices for the modified constraint bounds (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify c_U

Purpose

Modify upper bounds for nonlinear constraints. If idx is not given c U will be replaced.

Calling Syntax 

Prob = modify_c_U(Prob, c U, idx)

Description of Inputs

Prob Existing TOMLAB problem.
c_U New upper bounds for the nonlinear constraints.
idx Indices for the modified constraint bounds (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify x_0

Purpose

Modify starting point. If x_0 is outside the bounds an error will be returned. If idx is not given x 0 will be replaced.

Calling Syntax 

Prob = modify_x_0(Prob, x 0, idx)

Description of Inputs

Prob Existing TOMLAB problem.
x_0 New starting points.
idx Indices for the modified starting points (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify x_L

Purpose

Modify lower bounds for decision variables. If idx is not given x L will be replaced. x 0 will be shifted if needed.

Calling Syntax 

Prob = modify_x_L(Prob, x_L, idx)

Description of Inputs

Prob Existing TOMLAB problem.
x_L New lower bounds for the decision variables.
idx Indices for the modified lower bounds (optional).

Description of Outputs

Prob Modified TOMLAB problem.

modify x_U

Purpose

Modify upper bounds for decision variables. If idx is not given x U will be replaced. x 0 will be shifted if needed.

Calling Syntax 

Prob = modify_x_U(Prob, x_U, idx)

Description of Inputs

Prob Existing TOMLAB problem.
x_U New upper bounds for the decision variables.
idx Indices for the modified upper bounds (optional).

Description of Outputs

Prob Modified TOMLAB problem.

TomSym

For further information about TomSym, please visit http://tomsym.com/ - the pages contain detailed modeling examples and real life applications. All illustrated examples are available in the folder /tomlab/tomsym/examples/ in the TOMLAB installation. The modeling engine supports all problem types in TOMLAB with some minor exceptions.

A detailed function listing is available in TOMLAB Appendix C.

TomSym combines the best features of symbolic differentiation, i.e. produces source code with simplifications and optimizations, with the strong point of automatic differentiation where the result is a procedure, rather than an expression. Hence it does not grow exponentially in size for complex expressions.

Both forward and reverse modes are supported, however, reverse is default when computing the derivative with respect to more than one variable. The command derivative results in forward mode, and derivatives in reverse mode.

TomSym produces very efficient and fully vectorized code and is compatible with TOMLAB /MAD for situations where automatic differentiation may be the only option for parts of the model.

It should also be noted that TomSym automatically provides first and second order derivatives as well as problem sparsity patterns. With the use of TomSym the user no longer needs to code cumbersome derivative expressions and Jacobian/Hessian sparsity patterns for most optimization and optimal control problems.

The main features in TomSym can be summarized with the following list:

  • A full modeling environment in Matlab with support for most built-in mathematical operators.
  • Automatically generates derivatives as Matlab code.
  • A complete integration with PROPT (optimal control platform).
  • Interfaced and compatible with MAD, i.e. MAD can be used when symbolic modeling is not suitable.
  • Support for if, then, else statements.
  • Automated code simplification for generated models.
  • Ability to analyze most p-coded files (if code is vectorized).

Modeling

One of the main strength of TomSym is the ability to automatically and quickly compute symbolic derivatives of matrix expressions. The derivatives can then be converted into efficient Matlab code.

The matrix derivative of a matrix function is a fourth rank tensor - that is, a matrix each of whose entries is a matrix. Rather than using four-dimensional matrices to represent this, TomSym continues to work in two dimensions. This makes it possible to take advantage of the very efficient handling of sparse matrices in Matlab (not available for higher-dimensional matrices).

In order for the derivative to be two-dimensional, TomSym's derivative reduces its arguments to one-dimensional vectors before the derivative is computed. In the returned J , each row corresponds to an element of F , and each column corresponds to an element of X . As usual in Matlab, the elements of a matrix are taken in column-first order.

For vectors F and X , the resulting J is the well-known Jacobian matrix.

Observe that the TomSym notation is slightly different from commonly used mathematical notation. The notation used in tomSym was chosen to minimize the amount of element reordering needed to compute gradients for common expressions in optimization problems. It needs to be pointed out that this is different from the commonly used mathematical notation, where the tensor is flattened into a two-dimensional matrix as it is written (There are actually two variations of this in common use - the indexing of the elements of X may or may not be transposed).

For example, in common mathematical notation, the so-called self derivative matrix is a mn-by-mn square (or mm-by-nn rectangular in the non-transposed variation) matrix containing mn ones spread out in a random-looking manner. In tomSym notation, the self-derivative matrix is the mn-by-mn identity matrix.

The difference in notation only involves the ordering of the elements, and reordering the elements to a different notational convention should be trivial if tomSym is used to generate derivatives for applications other than for TOMLAB and PROPT.

Example: 

>> toms	y

>> toms 3x1 x

>> toms 2x3 A

>> f = (A*x).^(2*y)


f = tomSym(2x1): 

   (A*x).^(2*y)

>> derivative(f,A)

ans = tomSym(2x6):

   (2*y)*setdiag((A*x).^(2*y-1))*kron(x',eye(2))

In the above example, the 2x1 symbol f is differentiated with respect to the 2x3 symbol A. The result is a 2x6 matrix, representing .

The displayed text is not necessarily identical to the m-code that will be generated from an expression. For example, the identity matrix is generated using speye in m-code, but displayed as eye (Derivatives tend to involve many sparse matrices, which Matlab handles efficiently). The mcodestr command converts a tomSym object to a Matlab code string. 

>> mcodestr(ans)
ans =
setdiag((2*y)*(A*x).^(2*y-1))*kron(x',[1 0;0 1])

Observe that the command mcode and not mcodestr should be used when generating efficient production code.

Ezsolve

TomSym provides the function ezsolve, which needs minimal input to solve an optimization problem: only the objective function and constraints. For example, the miqpQG example from the tomlab quickguide can be reduced to the following: 

toms integer  x 
toms	      y

objective = -6*x  + 2*x^2  + 2*y^2  - 2*x*y;
constraints = {x+y<=1.9,x>=0, y>=0};

solution = ezsolve(objective,constraints)

Ezsolve calls tomDiagnose to determine the problem type, getSolver to find an appropriate solver, then sym2prob, tomRun and getSoluton in sequence to obtain the solution.

Advanced users might not use ezsolve, and instead call sym2prob and tomRun directly. This provides for the possibility to modify the Prob struct and set flags for the solver.

Usage

TomSym, unlike most other symbolic algebra packages, focuses on vectorized notation. This should be familiar to Matlab users, but is very different from many other programming languages. When computing the derivative of a vector-valued function with respect to a vector valued variable, tomSym attempts to give a derivative as vectorized Matlab code. However, this only works if the original expressions use vectorized notation. For example: 

toms 3x1 x
f = sum(exp(x));
g = derivative(f,x);

results in the efficient . In contrast, the mathematically equivalent but slower code: 

toms 3x1 x 
f = 0;
for i=1:length(x)
   f = f+exp(x(i));
end
g = derivative(f,x);

results in as each term is differentiated individually. Since tomSym has no way of knowing that all the terms have the same format, it has to compute the symbolic derivative for each one. In this example, with only three iterations, that is not really a problem, but large for-loops can easily result in symbolic calculations that require more time than the actual numeric solution of the underlying optimization problem.

It is thus recommended to avoid for-loops as far as possible when working with tomSym.

Because tomSym computes derivatives with respect to whole symbols, and not their individual elements, it is also a good idea not to group several variables into one vector, when they are mostly used individually. For example: 

toms 2x1 x
f = x(1)*sin(x(2));
g = derivative(f,x);

results in . Since x is never used as a 2x1 vector, it is better to use two independent 1x1 variables: 

toms a b
f = a*sin(b);
g = derivative(f,[a; b]);

which results in . The main benefit here is the increased readability of the auto-generated code, but there is also a slight performance increase (Should the vector x later be needed, it can of course easily be created using the code

Scaling variables

Because tomSym provides analytic derivatives (including second order derivatives) for the solvers, badly scaled problems will likely be less problematic from the start. To further improve the model, tomSym also makes it very easy to manually scale variables before they are presented to the solver. For example, assuming that an optimization problem involves the variable x which is of the order of magnitude 1e6, and the variable y, which is of the order of 1e - 6, the code: 

toms xScaled  yScaled 
x = 1e+6*xScaled;
y = 1e-6*yScaled;

will make it possible to work with the tomSym expressions x and y when coding the optimization problem, while the solver will solve for the symbols xScaled and yScaled, which will both be in the order of 1. It is even possible to provide starting guesses on x and y (in equation form), because tomSym will solve the linear equation to obtain starting guesses for the underlying xScaled and yScaled.

The solution structure returned by ezsolve will of course only contain xScaled and yScaled, but numeric values for x and y are easily obtained via, e.g. subs(x,solution).

SDP/LMI/BMI interface

An interface for bilinear semidefinite problems is included with tomSym. It is also possible to solve nonlinear problems involving semidefinite constraints, using any nonlinear solver (The square root of the semidefinte matrix is then introduced as an extra set of unknowns).

See the examples optimalFeedbackBMI and example sdp.

Interface to MAD and finite differences

If a user function is incompatible with tomSym, it can still be used in symbolic computations, by giving it a "wrapper". For example, if the cos function was not already overloaded by tomSym, it would still be possible to do the equivalent of cos(3*x) by writing feval('cos',3*x).

MAD then computes the derivatives when the Jacobian matrix of a wrapped function is needed. If MAD is unavailable, or unable to do the job, numeric differentiation is used.

Second order derivatives cannot be obtained in the current implementation.

It is also possible to force the use of automatic or numerical differentiation for any function used in the code. The follow examples shows a few of the options available: 

toms x1 x2 
alpha = 100;

%   1.  USE  MAD  FOR  ONE  FUNCTION.
%   Create  a wrapper  function. In  this case we  use sin, but  it could  be any
%  MAD   supported  function.
y = wrap(struct('fun','sin','n',1,'sz1',1,'sz2',1,'JFuns','MAD'),x1/x2);
f = alpha*(x2-x1^2)^2 + (1-x1)^2 + y;

%   Setup and solve  the  problem 
c = -x1^2 - x2;
con = {-1000  <= c <= 0
    -10  <= x1 <= 2
    -10  <= x2 <= 2};

x0 = {x1  == -1.2 
    x2 == 1};

solution1 = ezsolve(f,con,x0);

%   2.  USE  MAD  FOR  ALL FUNCTIONS. 
options = struct; 
options.derivatives = 'automatic';
f = alpha*(x2-x1^2)^2 + (1-x1)^2 + sin(x1/x2);
solution2 = ezsolve(f,con,x0,options);

%   3.  USE  FD  (Finite  Differences) FOR  ONE  FUNCTIONS.
%   Create  a new wrapper  function. In  this case we  use sin, but  it could  be
%   any function since  we  use numerical derivatives.
y = wrap(struct('fun','sin','n',1,'sz1',1,'sz2',1,'JFuns','FDJac'),x1/x2);
f = alpha*(x2-x1^2)^2 + (1-x1)^2 + y;
solution3 = ezsolve(f,con,x0);

%   4.  USE  FD  and MAD  FOR  ONE  FUNCTION    EACH.
y1 = 0.5*wrap(struct('fun','sin','n',1,'sz1',1,'sz2',1,'JFuns','MAD'),x1/x2);
y2 = 0.5*wrap(struct('fun','sin','n',1,'sz1',1,'sz2',1,'JFuns','FDJac'),x1/x2);
f = alpha*(x2-x1^2)^2 + (1-x1)^2 + y1 + y2;
solution4 = ezsolve(f,con,x0);

%   5.  USE  FD  FOR  ALL FUNCTIONS. 
options = struct; 
options.derivatives = 'numerical';
f = alpha*(x2-x1^2)^2 + (1-x1)^2 + sin(x1/x2);
solution5 = ezsolve(f,con,x0,options);

%   6.  USE MAD FOR OBJECTIVE,   FD FOR CONSTRAINTS
options = struct; 
options.derivatives = 'numerical'; 
options.use_H	= 0; 
options.use_d2c = 0;
options.type	= 'con';
Prob = sym2prob(f,con,x0,options);
madinitglobals; Prob.ADObj = 1; Prob.ConsDiff = 1;
Result  = tomRun('snopt', Prob,  1);
solution6 = getSolution(Result);

Simplifications

The code generation function detects sub-expressions that occur more than once, and optimizes by creating tem- porary variables for those since it is very common for a function to share expressions with its derivative, or for the derivative to contain repeated expressions.

Note that it is not necessary to complete code generation in order to evaluate a tomSym object numerically. The subs function can be used to replace symbols by their numeric values, resulting in an evaluation.

TomSym also automatically implements algebraic simplifications of expressions. Among them are:

  • Multiplication by 1 is eliminated: 1 * A = A
  • Addition/subtraction of 0 is eliminated: 0 + A = A
  • All-same matrices are reduced to scalars: [3; 3; 3] + x = 3 + x
  • Scalars are moved to the left in multiplications: A * y = y * A
  • Scalars are moved to the left in addition/subtraction: A - y = -y + A
  • Expressions involving element-wise operations are moved inside setdiag: setdiag(A)+setdiag(A) = setdiag(A+ A)
  • Inverse operations cancel: sqrt(x)2 = x
  • Multiplication by inverse cancels: A * inv(A) = eye(size(A))
  • Subtraction of self cancels: A - A = zeros(size(A))
  • Among others...

Except in the case of scalar-matrix operations, tomSym does not reorder multiplications or additions, which means that some expressions, like (A+B)-A will not be simplified (This might be implemented in a later version, but must be done carefully so that truncation errors are not introduced).

Simplifications are also applied when using subs. This makes it quick and easy to handle parameterized problems. For example, if an intricate optimization problem is to be solved for several values of a parameter a, then one might first create the symbolic functions and gradients using a symbolic a, and then substitute the different values, and generate m-code for each substituted function. If some case, like a = 0 results in entire sub-expressions being eliminated, then the m-code will be shorter in that case.

It is also possible to generate complete problems with constants as decision variables and then change the bounds for these variables to make them "real constants". The backside of this is that the problem will be slightly larger, but the problem only has to be generated once.

The following problem defines the variable alpha as a toms, then the bounds are adjusted for alpha to solve the problem for all alphas from 1 to 100. 

toms x1 x2

%   Define  alpha  as a toms although it is a constant toms alpha

%   Setup and solve  the  problem
f = alpha*(x2-x1^2)^2 + (1-x1)^2;
c = -x1^2  - x2;
con = {-1000  <= c <= 0
    -10  <= x1 <= 2
    -10  <= x2 <= 2};
x0 = {x1  == -1.2; x2 == 1}; 

Prob = sym2prob(f,con,x0);

%  Now   solve  for alpha  = 1:100, while reusing x_0 
obj  = zeros(100,1);

for  i=1:100
    Prob.x_L(Prob.tomSym.idx.alpha) = i; 
    Prob.x_U(Prob.tomSym.idx.alpha) = i; 
    Prob.x_0(Prob.tomSym.idx.alpha) = i; 
    Result  = tomRun('snopt', Prob,  1); 
    Prob.x_0  = Result.x_k;
    obj(i) = Result.f_k;
end

Special functions

TomSym adds some functions that duplicates the functionality of Matlab, but that are more suitable for symbolic treatment. For example:

  • setDiag and getDiag - Replaces some uses of Matlab's diag function, but clarifies whether diag(x) means "create a matrix where the diagonal elements are the elements of x" or "extract the main diagonal from the matrix x".
  • subsVec applies an expression to a list of values. The same effect can be achieved with a for-loop, but subsVec gives more efficient derivatives.
  • ifThenElse - A replacement for the if ... then ... else constructs (See below).

If ... then ... else:

A common reason that it is difficult to implement a function in tomSym is that it contains code like the following: 

if x<2
    y = 0;
else
    y = x-2;
end

Because x is a symbolic object, the expression x < 2 does not evaluate to true or false, but to another symbolic object.

In tomSym, one should instead write: 

y = ifThenElse(x<2,0,x-2)

This will result in a symbolic object that contains information about both the "true" and the "false" scenario. However, taking the derivative of this function will result in a warning, because the derivative is not well-defined at x = 2.

The "smoothed" form: 

y = ifThenElse(x<2,0,x-2,0.1)

yields a function that is essentially the same for abs(x - 2) > 3 * 0.1, but which follows a smooth curve near x = 2, ensuring that derivatives of all orders exist. However, this introduces a local minimum which did not exist in the original function, and invalidates the convexity.

It is recommended that the smooth form ifThenElse be used for nonlinear problems whenever it replaces a dis- continuous function. However, for convex functions (like the one above) it is usually better to use code that helps tomSym know that convexity is preserved. For example, instead of the above if ThenElse(x < 2, 0, x - 2, 0.1), the equivalent max(0, x - 2) is preferred.

Procedure vs parse-tree

TomSym works with procedures. This makes it different from many symbolic algebra packages, that mainly work with parse-trees.

In optimization, it is not uncommon for objectives and constraints to be defined using procedures that involve loops. TomSym is built to handle these efficiently. If a function is defined using many intermediate steps, then tomSym will keep track of those steps in an optimized procedure description. For example, consider the code: 

toms x
y = x*x;
z = sin(y)+cos(y);

In the tomSym object z, there is now a procedure, which looks something like: 

temp = x*x;
result = cos(temp)+sin(temp);

Note: It is not necessary to use the intermediate symbol y. TomSym, automatically detects duplicated expressions, so the code would result in the same optimized procedure for z.

On the other hand, the same corresponding code using the symbolic toolbox: 

syms x
y = x*x;
z = sin(y)+cos(y);

results in an object z that contains , resulting in a double evaluation of .

This may seem like a small difference in this simplified example, but in real-life applications, the difference can be significant.

Numeric stability:

For example, consider the following code, which computes the Legendre polynomials up to the 100th order in tomSym (The calculation takes about two seconds on a modern computer). 

toms x 
p{1}=1; 
p{2}=x;
for i=1:99
   p{i+2} = ((2*i+1)*x.*p{i+1}-i*p{i})./(i+1);
end

Replacing "toms" by "syms" on the first line should cause the same polynomials to be computed using Mathwork's Symbolic Toolbox. But after a few minutes, when only about 30 polynomials have been computed, the program crashes as it fails to allocate more memory. This is because the expression grows exponentially in size. To circumvent the problem, the expression must be algebraically simplified after each step. The following code succeeds in computing the 100 polynomials using the symbolic toolbox. 

syms x 
p{1}=1; 
p{2}=x;
for i=1:99
   p{i+2} = simplify(((2*i+1)*x.*p{i+1}-i*p{i})./(i+1));
end

However, the simplification changes the way in which the polynomial is computed. This is clearly illustrated if we insert x = 1 into the 100th order polynomial. This is accomplished by using the command subs(p101,x,1) for both the tomSym and the Symbolic Toolbox expressions. TomSym returns the result 1.0000, which is correct. The symbolic toolbox, on the other hand, returns 2.6759e + 020, which is off by 20 orders of magnitude. The reason is that the "simplified" symbolic expressions involves subtracting very large numbers. Note: It is of course possible to get correct results from the Symbolic Toolbox using exact arithmetic instead of machine-precision floating-point, but at the cost of much slower evaluation.

In tomSym, there are also simplifications, for example identifying identical sub-trees, or multiplication by zero, but the simplifications are not as extensive, and care is taken to avoid simplifications that can lead to truncation errors. Thus, an expression computed using tomSym should be exactly as stable (or unstable) as the algorithm used to generate it.

Another example:

The following code, iteratively defines q as a function of the tomSym symbol x, and computes its derivative: 

toms x 
q=x;
for i=1:4
   q = x*cos(q+2)*cos(q);
end 
derivative(q,x)

This yields the following tomSym procedure: 

tempC3	= x+2;
tempC4	= cos(tempC3); 
tempC5	= x*tempC4; 
tempC10	= cos(x);
tempC12	= tempC10*(tempC4-x*sin(tempC3))-tempC5*sin(x);
tempC13	= tempC5*tempC10; 
tempC16	= tempC13+2; 
tempC17	= cos(tempC16); 
tempC18	= x*tempC17; 
tempC24	= cos(tempC13);
tempC26	= tempC24*(tempC17-x*(sin(tempC16)*tempC12))-tempC18*(sin(tempC13)*tempC12);
tempC27	= tempC18*tempC24; 
tempC30	= tempC27+2; 
tempC31	= cos(tempC30); 
tempC32	= x*tempC31; 
tempC38	= cos(tempC27);
tempC40	= tempC38*(tempC31-x*(sin(tempC30)*tempC26))-tempC32*(sin(tempC27)*tempC26);
tempC41	= tempC32*tempC38; 
tempC44	= tempC41+2; 
tempC45	= cos(tempC44);
out 	= cos(tempC41)*(tempC45-x*(sin(tempC44)*tempC40))-(x*tempC45)*(sin(tempC41)*tempC40);

Now, complete the same task using the symbolic toolbox: 

syms x q=x;
for i=1:4
   q = x*cos(q+2)*cos(q);
end 
diff(q,x)

This yields the following symbolic expression: 

cos(x*cos(x*cos(cos(x+2)*x*cos(x)+2)*cos(cos(x+2)*x*cos(x))+2)*cos(x*cos(cos(x+2)*x*cos(x)+2)*... 

cos(cos(x+2)*x*cos(x)))+2)*cos(x*cos(x*cos(cos(x+2)*x*cos(x)+2)*cos(cos(x+2)*x*cos(x))+2)*cos(x*... 

cos(cos(x+2)*x*cos(x)+2)*cos(cos(x+2)*x*cos(x))))-x*sin(x*cos(x*cos(cos(x+2)*x*cos(x)+2)*cos(...

...

and 23 more lines of  code.

And this example only had four iterations of the loop. Increasing the number of iterations, the Symbolic toolbox expressions quickly becomes unmanageable, while the tomSym procedure only grows linearly.

Problems and error messages

  • Warning: Directory c:\Temp\tp563415 could not be removed (or similar). When tomSym is used to automatically create m-code it places the code in a temporary directory given by Matlab's tempname function. Sometimes Matlab chooses a name that already exists, which results in this error message (The temporary directory is cleared of old files regularly by most modern operating systems. Otherwise the temporary Matlab files can easily be removed manually).
  • Attempting to call SCRIPT as a function (or similar). Due to a bug in the Matlab syntax, the parser cannot know if f (x) is a function call or the x:th element of the vector f. Hence, it has to guess. The Matlab parser does not understand that toms creates variables, so it will get confused if one of the names is previously used by a function or script (For example, "cs" is a script in the systems identification toolbox). Declaring toms cs and then indexing cs(1) will work at the Matlab prompt, but not in a script. The bug can be circumvented by assigning something to each variable before calling toms.

Example

A TomSym model is to a great extent independent upon the problem type, i.e. a linear, nonlinear or mixed-integer nonlinear model would be modeled with about the same commands. The following example illustrates how to construct and solve a MINLP problem using TomSym. 

Name='minlp1Demo - Kocis/Grossman.';

toms 2x1 x
toms 3x1 integer  y

objective = [2  3 1.5  2 -0.5]*[x;y];

constraints = { ... 
   x(1) >= 0,  ... 
   x(2) >= 1e-8,   ... 
   x <= 1e8,  ...
   0 <= y <=1,  ...
   [1 0 1 0 0]*[x;y] <= 1.6, ...
   1.333*x(2) + y(2) <= 3,  ... 
   [-1 -1  1]*y <= 0,  ...
   x(1)^2+y(1) == 1.25, ... 
   sqrt(x(2)^3)+1.5*y(2) == 3,  ...
};

guess = struct('x',ones(size(x)),'y',ones(size(y)));
options = struct;
options.name  = Name;
Prob = sym2prob('minlp',objective,constraints,guess,options); 
Prob.DUNDEE.optPar(20) = 1;
Result  = tomRun('minlpBB',Prob,2);

The TomSym engine automatically completes the separation of simple bounds, linear and nonlinear constraints.

Retrieved from "http://tomwiki.com/index.php?title=TOMLAB_Defining_Problems_in_TOMLAB&oldid=2435"