TomSym
TOMSYM is a symbolic manipulation package which is part of the TOMLAB optimization software suite. It simplifies the work of defining optimization problems by automatically generating derivatives, m-code and other things needed for the numeric solvers to work as efficiently as possible.
TOMSYM is also the foundation on which PROPT is built. We recommend that PROPT users read through this guide first, before continuing on to the PROPT manual.
Solving optimization problems
An optimization problem is defined by an objective function and a set of constraints. Problems are usually classified depending on the properties of the objective and constraints, and different types of numeric techniques are used for different types of problems. For example, solvers usually handle quadratic objectives and linear constraints more efficiently than general nonlinear ones.
TOMSYM provides the ezsolve function for solving any type of optimization problem. It examines the objective and constraints to determine the problem type, and then automatically selects the numeric solver which it deems most apropriate.
For example, here we find a minimum to a second degree polynomial in two variables
>> toms a b
>> ezsolve(a^2 + 4*b^2 + 3*a - 2*b - a*b)
Problem type appears to be: qp
Time for symbolic processing: 0.040062 seconds
Starting numeric solver
=====* * * ===================================================================
TOMLAB - Tomlab Optimization Inc. Development license 999001. Valid to 2013-02
===============================================================================
Problem: 1: f_k -2.266666666666666607
f(x_0) 0.000000000000000000
Solver: CPLEX. EXIT=0. INFORM=1.
CPLEX Barrier QP solver
Optimal solution found
FuncEv 2 GradEv 2 ConstrEv 2 Iter 2
CPU time: 0.010000 sec. Elapsed time: 0.002605 sec.
ans =
a: -1.4667
b: 0.0667
Looking at the output from the solver, we see that ezsolve determined the problem type to be QP (quadratic programming) and selected CPLEX at the solver. The problem was then solved numerically and the solution is returned as a matlab structure with fields corresponding to the symbol names.
Most solvers handle the linear constraints separately from the nonlinear ones. With TOMSYM, this happens transparently. The linear constraints are collected and passed to the solver in the form of a (sparse) matrix, while the nonlinear ones are passed as m-code, together with m-code for their derivatives.
Symbolic matrices
With TOMSYM, a single symbol can be a scalar, vector or matrix. Expressions that involve matrix symbols are stored only once. This behaviour is defferent from many symbolic algebra packages, where symbolic arrays are treated as arrays of symbols (i.e. one symbolic expression is stored for each position in the array).
Symbolic matrices can be manipulated just like normal MATLAB matrices. For example, * gives the matrix product, and .* gives the element-wise product of two matrices.
>> A = tom('A',2,2)
A = tomSym(2x2):
A
>> b = tom('b',2,1)
b = tomSym(2x1):
b
>> y = A*b
y = tomSym(2x1):
A*b
In this example A is 2-by-2 matrix symbol, b is a 2-by-1 column vector symbol, and the matrix expression y is also a 2-by-1 column vector.
Taking the derivative of a vector expression with respect to a vector symbol gives the so-called Jacobian matrix.
>> derivative(y,b) ans = tomSym(2x2): A
Derivatives involving matrices are computed as if the matrices were vectors, with elements taken in column-first order.
>> derivative(y,A) ans = tomSym(2x4): kron(b',eye(2))
A consequence of this definition is that the derivative of a matrix with respect to itself is the identity matrix.
>> full(derivative(A,A))
ans =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
Vectorized code
In languages such as C or FORTRAN, "loops" are very common occurences. Matlab, on the other hand, uses a syntax where loops can usually be replaced by vectorized expressions. An example is the matrix product, which can be coded in C using three nested for loops, but in Matlab is accomplished by simply using the operator *.
tomSym is intended for use with vectorized code. It supports all the functions that are typically needed to generated vectorized expressions, such as sum, diff, repmat, sparse, etc.
Complex numbers
The solvers in the TOMLAB suite work with real numbers. Still, some equations are more easily expressed using complex numbers.
The tomCmplx class makes it possible to represent symbolic complex expressions as pairs of ordinary tomSym expressions. For example:
>> toms a b >> z = tomCmplx(a,b) z = tomCmplx(1x1): Real part: a Imaginary part: b >> z2 = z.*z z2 = tomCmplx(1x1): Real part: a*a-b*b Imaginary part: a*b+b*a
The set of operations allowed for tomCmplx expressions is currently very limited. Element-wise addition, subtraction, multiplication and division are supported, as are real(), imag(), abs() and angle() which return real-valued tomSym expressions.
If a tomCmplx is used in an equation, then a cell array of two equations is generated; One for the real and one for the imaginary component. A complex number can never be used in an inequality, so <= and >= are not supported for tomCmplx.
What tomSym cannot do
There are a few things that are possible in Matlab code, but which tomSym currently cannot do. Below we list the most common issues, and their work-arounds.
Variable-size matrixes
The size of a tomSym object cannot depend on the value of another tomSym object. This has to do with the fact that numeric solvers typically allocate memory only at the start of an optimization run, and then expect the number of unknowns to remain constant. (This may change in the future.)
Indexes into double arrays
Replacing values in arrays
Conditional execution
N-dimensional arrays for N > 2
TOMSYM works with scalars, vectors, and two-dimensional matrices. Although Matlab can handle arrays of any dimension, TOMSYM cannot. The work-around is to use tomArray, which transforms tomSym symbols into arrays of any dimension. (TomArray also allows for a GAMS-style indexing which is a bit different from Matlab syntax, but convenient when formulating optimization problems involving high-dimension arrays.)
Note that TOMSYM supports interpn with high-dimensional numeric data, although the symbolic inputs must be of dimension less than two.
Limits
TOMSYM evaluates expressions by generating matlab code. (It does some mathematical simplification to improve efficiency, but only to a limited extent.) If an expression evaluates to 0/0 or 0*Inf, then NaN (Not a Number) will be returned. When the solver encounters a NaN, it will often just quit (although some solvers include back-tracking algorithms.)
This means that expressions such as a*log(a) can not be used if a might become zero. In this case it might be prudent to add a constraint a >= 1e-8 (or similar) to help the solver stay out of trouble.
Interfaces
TOMSYM is most commonly used for setting up optimization problems for the TOMLAB suite of numeric solvers. However, it is also possible to use it by itself to generate m-code or to interface it with your own code.
Generating m-code with TOMSYM
It is possible to convert any TOMSYM object into MATLAB code via the mcode and mfile commands. This makes it possible, for example, to use derivatives computed by TOMSYM in your own applications outside of TOMLAB.
Because it is very inefficient to store large matrices as code, the mcode/mfile commands also return a cell array named tempD which contains the numeric data from the TOMSYM object. The tempD array must be provided whenever the code is executed.
Connection to TOMLAB
When TOMLAB is used without TOMSYM, optimization problems are represented via a Prob structure, which typically contain function handles to the objective and constraint functions, and ideally also their derivatives. TOMSYM uses a symbolic representation of the objective and constraints to create this Prob structure, including the m-code needed, and pass it directly to the solver. This and enables the user to focus on the modelling task rather than tedious implementation details such as coding up the derivative for each nonlinear constraint.
Advanced users may extract the Prob structure from TOMSYM using the sym2prob function, and then manipulate this structure further before calling the solver. This makes it possible, for example, to optimize the autogenerated code, or set special solver flags.
When the same problem needs to be re-solved many times for varying values of a parameter, the fastest way is usually to create the Prob structure once, and then change the parameter value using setparameter
Connection to PROPT
PROPT is a solver for dynamic optimization problems which uses a collocation mehtod. PROPT is based on TOMSYM, and uses the same symbolic engine. The output of PROPT's collocate family of functions is ordinary TOMSYM arrays which means that PROPT and TOMSYM code can be mixed freely.
Improving speed
Speed is often important when solving optimization problems. The main advantage of tomSym is that it allows for high overall speed. The time it takes from the moment you start coding until the solution presented should be as short as possible. Often, much more time is spent coding than numerically solving the problem, so emphasis should be placed on writing code that is easy to read and debug. Sometimes, however, the time used by the numerical solver is critical, and this section gives som hints on how to make it run faster.
Pre-process the symbolic problem
See tomsym_mpc.m
Avoid loops
Manipulating tomSym expressons inside a loop can be very slow. For example, the following loop generates a long symbolic expression.
>> v = tom('v',1,10);
>> s = 0;
>> for i = 1:length(v); s = s+v(i); end
>> disp(s);
((((((((v(1)+v(2))+v(3))+v(4))+v(5))+v(6))+v(7))+v(8))+v(9))+v(10)
In this case, the loop can simply be replaced by s = sum(v), which will run much faster.
Use separate symbols for separate things
It is common in some optimization frameworks to use a single vector x for all the optimization variables, and then using indexes to represent different things. (x(1) is position in m, x(2) is speed in m/s, etc.)
Doing this in tomSym is a bad idea for several reasons. The generated code will be less efficient, because tomSym tries to generated derivatives with respect to the whole vector but it is only used one element at a time. Automatic scaling will not work as expected, becuase tomSym tries to scale the entire vector by one factor, when the individual elements may have different magnintude.
Also, it is much easier to read code that uses x for position and v for speed, than code that uses x for both.
Don't vectorize expressions that don't belong together
It is faster say { x >= 0, cos(x) >= 0.5 } than [x; cos(x)] >= [ 0; 0.5 ].
In the first case we have one linear constraint and one nonlinar constraint. The solver will handle these separately, and linear constraints are treated much more efficiently than nonlinear ones. In the second case we have a nonlinear constraint with two components. TomSym is unable to separate the linear component, so both are handled by the nonlinear part of the solver.
Avoid using abs and max
It is much more efficient to say { v >= -1, v <= 1 } than { max(abs(v)) <= 1 }. Again, this is because linear constraints are treated more efficiently than nonlinear ones. For the same reason it is faster to minimize y under the constraint { v >= -y, v <= y } than to minimize max(abs(x)).
The rewriteV function (called by ezsolve) attempts to get rid of instances of abs and max, but it is only able to handle simple cases.
Don't put equations in the objective
Using an objective that requires solving equations, such as sum((Ab).^2)) is very inefficient. A better way to handle this is to rewrite the objective as sum(x.^2) and add the constraint A*x==b. This formulation gives much more information to the solver, and avoids having to compute the derivative of the solution to an equation (which is numerically unstable.)
For this reason, functions like fzero should never be used with tomSym. The function tomfzero can often be used in its place.
Balance the scale
Equations like a/1000000 + b ==1 don't work well with numeric solvers, because a and b will have very different magnitudes in the solution.
One way to overcome this problem is to set OPTIONS.scale = 'auto', which will cause ezsolve to try to guess a scaling factor for each unknown. This works best if each variable has been given reasonable upper and lower bounds.
Another way to fix the problem is to use scaled variable from the start. For example, if we use the scaled variable as, and let a equal 1000000*as, then the equation becomes as + b ==1, which is well-scaled.
Multiply rather than divide
Multiplication is "less nonlinear" than division. Hence, equations like m*x ==f typically result in better convergence than x ==f/m. This is particularly true for matrices: M*x ==f can be much faster than x ==Mf.
This method also prevents the problem where the solver shuts down because of a division-by-zero. (However, if m is zero then f== is a solution, so watch out!)
Supported functions
Below is a list of functions that are overloaded for tomSym. That is: these functions accept symbolic input and then return symbolic expressions whose derivatives can be computed symbolically.
Most of these are standard Matlab functions, but some are special functions that are mostly used internally by tomSym for efficiency reasons.
>> help tomSym/Contents
@TOMSYM
Files
abs - tomSym/abs - Overloaded function
acos - tomSym/acos - Overloaded function
acosh - tomSym/acosh - Overloaded function
acot - tomSym/acot - Overloaded function
acoth - tomSym/acoth - Overloaded function
acsc - tomSym/acsc - Overloaded function
acsch - tomSym/acsch - Overloaded function
angle - tomSym/angle - Overloaded function
asec - tomSym/asec - Overloaded function
asech - tomSym/asech - Overloaded function
asin - tomSym/asin - Overloaded function
asinh - tomSym/asinh - Overloaded function
atan - tomSym/atan - Overloaded function
atan2 - tomSym/atan2 - Overloaded function
atanh - tomSym/atanh - Overloaded function
bsxfun - tomSym/bsxfun - Overloaded function
cat - tomSym/cat - Overloaded function
ceil - tomSym/ceil - Overloaded function
char - tomSym/char - Convert tomSym to char.
complementary - complementary - overloaded function
conj - tomSym/conj - Overloaded function
constpart - tomSym/constpart - The zeroth term of the McLau
cos - tomSym/cos - Overloaded function
cosh - tomSym/cosh - Overloaded function
cot - tomSym/cot - Overloaded function
csc - tomSym/csc - Overloaded function
ctranspose - tomSym/ctranspose - Overloaded operator
cumsum - tomSym/cumsum - Overloaded function
dblquad - tomSym/dblquad - Numeric double quadrature
derivative - tomSym/derivative - The symbolic derivative of
derivatives - tomSym/derivatives - The symbolic derivatives o
det - tomSym/det - Overloaded function
diag - tomSym/diag - Overloaded function
diff - tomSym/diff - Difference - Overloaded function
DiracDelta - tomSym/DiracDelta - Placeholder derivative of a
DiracDeltas - tomSym/DiracDeltas - Placeholder derivative of
disp - tomSym/disp - Command window display of a tomSy
display - tomSym/display - Command window display of a to
double - tomSym/double - Convert tomSym constant to doub
end - tomSym/end - Overloaded - Get the last index of
eps - tomSym/rem - Overload the remainder function
eq - tomSym/eq - Overloaded ==operator
erf - tomSym/erf - Overloaded function
erfc - tomSym/erfc - Overloaded function
erfcinv - tomSym/erfcinv - Overloaded function
erfinv - tomSym/erfinv - Overloaded function
estpattern - tomSym/estpattern - Estimate the sparsity patte
eval - tomSym/eval - Evaluate a tomSym object in the c
exp - tomSym/exp - Overloaded function
extractConstraints - extrectContraitns - Move constraints from subje
eye - tomSym/eye - Overloaded function
ezplot - tomSym/ezplot - Function plot
feval - tomSym/feval - tomSym-compatible function call.
fix - tomSym/fix - Overloaded function
fliplr - tomSym/fliplr - Overloaded function
flipud - tomSym/flipud - Overloaded function
floor - tomSym/floor - Overloaded function
full - tomSym/full - Overloaded function
fzero - tomSym/fzero - same as tomfzero
gamma - tomSym/gamma - Overloaded function
gammainc - tomSym/gammainc - Overloaded function
gammaln - tomSym/gammaln - Overloaded function
ge - tomSym/ge - Overloaded >= operator
getdiag - tomSym/getdiag - Getdiag for tomSym.
gt - tomSym/gt - Overloaded > operator
holdInterpolationMatrix - holdInterpolationMatrix - Value matrix for inte
horzcat - tomSym/horzcat - Overload the [] operator
hownonlinear - hownonlinear - overloaded function
ifThenElse - tomSym/ifThenElse - Symbolic if/then/else
imag - tomSym/imag - Overloaded function
interp1 - tomSym/interp1 - Overloaded function
interp1h - tomSym/interp1h - Overloaded function
interp1l - tomSym/interp1l - Overloaded function
interp1s - tomSym/interp1s - Overloaded function
interp1sDot - tomSym/interp1sDot - Derivative of interp1s, ov
interp2 - TOMSYM/INTERP2 - overloaded function
interpn - TOMSYM/INTERPN - overloaded function
inv - tomSym/inv - Overloaded function
isdependent - tomSym/isdependent - Determine if f is dependen
isempty - tomSym/isempty - Overloaded function
isequal - tomSym/isequal - Overloaded function
isfinite - tomSym/isfinite - Overloaded function
isinf - tomSym/isinf - Overloaded function
isnan - tomSym/isnan - Overloaded function
isreal - tomSym/isreal - Overloaded function
istomsymbol - tomSym/istomsymbol - Returns "true" if a is a t
kkt1 - kkt1 - First order Karush-Kuhn-Tucker condition
kron - tomSym/kron - Overload the Kronecker product
ldivide - tomSym/ldivide - Overload the .\ left division
le - tomSym/ge - Overloaded <= operator
length - tomSym/length - Get the length of a tomSym obje
linearInterpolationMatrix - linearInterpolationMatrix - Value matrix for in
linspace - tomSym/linspace - Overloaded function
log - tomSym/log - Overloaded function
log10 - tomSym/log10 - Overloaded function
log2 - tomSym/log2 - Overloaded function
logm - tomSym/logm - Overloaded function
lookup - tomSym/lookup - Overloaded function
lt - tomSym/lt - Overloaded < operator
mat2str - tomSym/mat2str - Convert tomSym to a string.
max - tomSym/max - Overloaded
maxIndicator - maxIndicator - Overloaded function
mcode - tomSym/mcode - Generate m-code from a tomSym ob
mcodestr - tomSym/mcodestr - Convert tomSym to m-code
MI - tomSym/MI - Matrix Inequality
min - tomSym/min - Overloaded
minus - tomSym/minus - Overload the minus operator
mldivide - tomSym/mrdivide - Overload the left divide oper
mod - tomSym/mod - Overloaded function
monofun - monofun - Overloaded function
monoinv - monoinv - Overloaded function
mpower - tomSym/mpower - Overload the mpower operator
mrdivide - tomSym/mrdivide - Overload the divide operator
mtimes - tomSym/mtimes - Overload the multiplication ope
nnz - tomSym/spy - Number of nonzero matrix elements
nOperands - tomSym/nOperands - Get number of operands from
norm - tomSym/norm - Overloaded function
not - tomSym/not - Overloaded function
num2str - tomSym/num2str - Convert tomSym to a string.
numel - tomSym/numel - Returns prod(size(a)) for tomSym
operand - tomSym/operand - Get an operand from a tomSym
operands - tomSym/operands - Get all operands from a tomSy
operator - tomSym/operator - Get the operator from a tomSy
pattern - tomSym/pattern - The sparsity pattern of a tomS
plus - tomSym/plus - Overload the plus operator
positiveSemidefinite - tomSym/positiveSemidefinite - Overloaded functi
power - tomSym/power - Overload the mpower operator
ppnval - TOMSYM/PPNVAL - Overloaded function
ppval - tomSym/ppval - Overloaded function
prod - tomSym/prod - Overloaded function
prodJ1 - tomSym/prodJ1 - overloaded function
prodJ1J1 - tomSym/prodJ1J1 - overloaded function
psi - tomSym/psi - Overloaded function
quad - tomSym/quad - Numeric quadrature
rdivide - tomSym/rdivide - Overload the ./ division opera
real - tomSym/real - Overloaded function
rem - tomSym/rem - Overload the remainder function
repmat - tomSym/repmat - Overloaded function
reshape - tomSym/reshape - Overloaded function
rewriteV - rewriteV - Rewrite an optimization problem to a
round - tomSym/round - Overloaded function
scalecolumns - tomSym/scalecolumns - Overloaded function
scalerows - tomSym/scalerows - Overloaded function
sec - tomSym/sec - Overloaded function
setdiag - tomSym/setdiag - Setdiag for tomSym.
setSymmetric - tomSym/setSymmetric - Create a symmetric matrix
sign - tomSym/sign - Overloaded function
sin - tomSym/sin - Overloaded function
sinh - tomSym/sinh - Overloaded function
size - tomSym/size - Get the size of a tomSym object
smplus - tomSym/smplus - Scalar + matrix addition
smtimes - tomSym/smtimes - Scalar x matrix multiplication
sparse - tomSym/sparse - Overloaded function
spdiags - tomSym/spdiags - Overloaded function
spower - tomSym/spower - Element-wise power to a scalar
spy - tomSym/spy - Visualize sparsity pattern for a t
sqrt - tomSym/sqrt - Overloaded function
sub2ind - tomSym/sub2ind - Overloaded function.
submatrix - tomSym/submatrix - Overloaded function
subsasgn - tomSym/subsasgn - Subscripted assignemnt, overl
subsindex - tomSym/subsindex - Not possible
subsref - tomSym/subsref - Object index lookup
subststruct - tomSym/subststruct - Substitue tomSym symbols f
subsymb - tomSym/subsymb - Get a sub-symbol from a tomSym
sum - tomSym/sum - Overloaded function
sym - tomSym/sym - convert a tomSym object to a symbo
sym2eig - sym2eig Convert symbolic constraints into an ei
sym2prob - sym2prob - Compile symbolic function/constraint
symbols - tomSym/symbols - List all symbols used in a tom
tan - tomSym/tan - Overloaded function
tanh - tomSym/tanh - Overloaded function
times - tomSym/times - Overload the .* multiplication o
tomCmp - tomCmp - Test if the operator of a tomSym is of
tomnumeric - tomSym/tomnumeric - isnumeric for tomSym
tomSym - tomSym/tomSym - Class constructor
tomUnGlobalize - tomSym/tomUnGlobalize - undo the effect of tomG
trace - tomSym/trace - Overloaded function
transpose - tomSym/transpose - Overloaded operator
tril - tomSym/tril - Overloaded function
triu - tomSym/triu - Overloaded function
uminus - tomSym/uminus - Overloaded function
uplus - tomSym/uplus - Overloaded function
vec - tomSym/vec - Transform a tomSym object into a c
vertcat - tomSym/vertcat - Overload the [;] operator
wrap - tomSym/wrap - Overloaded
wrapJ - tomSym/wrapJ - Overloaded
In addition to the functions on this list, there are many functions that are supported because they consist entirely of other functions that are supported. For example, it is possible that functions which you have written yourself already work with tomSym. The easiest way to find out if a function is supported is to try to use it. If you don't get an error message, then it works!
We constantly add new functions to TOMSYM. If there's a function that you think would be useful, let us know and we might include it in a future TOMLAB release.
