TomSym: Difference between revisions
(→Frequently Asked Questions: moved to TomSym FAQ) |
|||
Line 141: | Line 141: | ||
The set of operations allowed for tomCmplx expressions is currently very limited. Element-wise addition, subtraction, multiplication and | The set of operations allowed for tomCmplx expressions is currently very limited. Element-wise addition, subtraction, multiplication, division and exponentiation are supported, as are <tt>real()</tt>, <tt>imag()</tt>, <tt>abs()</tt> and <tt>angle()</tt> 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 <tt><=</tt> and <tt>>=</tt> are not supported for tomCmplx. | 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 <tt><=</tt> and <tt>>=</tt> are not supported for tomCmplx. (But inequalities can be used on the real and imaginary parts, as those are normal tomSym expressions.) | ||
==What tomSym cannot do== | ==What tomSym cannot do== |
Revision as of 10:01, 26 July 2011
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, division and exponentiation 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. (But inequalities can be used on the real and imaginary parts, as those are normal tomSym expressions.)
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.