# TomSym compared to the Symbolic Toolbox

We look at the optimization problem of maximizing the revenue of a hydroelectric dam. This problem was published by By Seth DeLand of the MathWorks, as an example of how to combine their optimization and symbolic toolboxes. [1] We solve the same probem using tomSym and the TOMLAB solvers, noting a significant speed difference, with tomSym running thousands of times faster than the Symbolic Toolbox.

## Problem description

The dynamics of the dam is described by the following equations.

Electricity(t) = TurbineFlow(t-1)*[½ k1(Storage(t) + Storage(t-1)) + k2]

Storage(t) = Storage(t-1) + ∆t * [inFlow(t-1) - spillFlow(t-1) - turbineFlow(t-1)]

The problem is to maximize revenue (electricity times price), subject to the following constraints:

- 0 < Turbine Flow < 25,000 CFS
- 0 < Spill Flow
- Max change in Turbine and Spill Flow < 500 CFS
- Combined Turbine and Spill Flow > 500 CFS
- 50000 < Reservoir Storage < 100000 Acre-Feet
- The final reservoir level must equal the initial level (90000 AF)

The code and data used by DeLand can be downloaded from Matlab Central.[2]

## tomSym code

To avoid attaching data files to our example, we will use an approximation to this data, as generated by the following Matlab code:

```
N = 480;
t = (0:N-1)';
inFlow = reshape(repmat(10*(107+[0 0 0 -3 3 2 2 -8 9 5 0 1 -6 -4 5 0 -6 -6 -17 -7]),24,1),N,1);
price = 46 + t./57 - 3*cos((t+46)/43) - 4*sin(pi/12*(t+4)) - 4*sin(pi/6*(t+1));
```

A few other consants are defined by DeLand as follows:

```
stor0 = 90000; % initial vol. of water stored in the reservoir (Acre-Feet)
k1 = 0.00003; % K-factor coefficient
k2 = 9; % K-factor offset
MW2kW = 1000; % MW to kW
C2A = 1.98347/24; % Convert from CFS to AF/HR
```

One the input data is defined, this is the entire code used to set up and solve the problem with TOMLAB:

```
turbineFlow = tom('turbineFlow',N,1);
spillFlow = tom('spillFlow',N,1);
outFlow = turbineFlow + spillFlow;
storage = stor0 + C2A*cumsum(inFlow - outFlow);
electricity = turbineFlow.*(0.5*k1*(storage + [stor0; storage(1:end-1)]) + k2)./MW2kW;
revenue = price'*electricity;
objective = -revenue; % Minimization objective
constraints = {
0 <= turbineFlow <= 25000
0 <= spillFlow
-500 <= diff(outFlow) <= 500
outFlow >= 500
50000 <= storage <= 100000
storage(end) == stor0
};
guess = [];
options = struct;
options.solver = 'knitro'; % also try: 'snopt'
options.Prob.KNITRO.options.ALG=3;
options.scale = 'auto'
solution = ezsolve(objective, constraints, guess, options);
```

## Solution

Although the Hessian is constant, this is *not* a quadratic programming.
The Hessian is not positive semidefinite, and the objective is therefore not convex. There may exist
local optima that are not as good as the global one. For this reason, the QP solvers in the
TOMLAB suite return error messages, and a nonlinear solver had to be used instead.

The solution for the approximate input data is very similar to the one presented by DeLand. With identical input data, the solutions are of course nearly identical (within the tolerances of the solvers).

Interestingly, the objective value reported by the quadprog solver does not match the value obtained from substituting the solution back into the objective function. The command

`(subs(Rtot,X,x3)-fval3)./fval3`

computes the relative error. This error should be zero, or on the order of machine precision: 1e-16, but it is on the order of 1e-4. This seems to be due to rounding errors in the gradient and Hessian as computed by the Symbolic Toolbox. The value obtained by "subs(Rtot,X,x3)" is very close to the optimum value returned by the TOMLAB solvers (within 1e-10 or so).

## Code comparison

One thing that is immediately noticeable when comparing the above code to DeLand's is that it is much shorter. After removing comments, blank lines and input data generation, the above code is about half the length of the part of DeLand's code that is needed for setting up and solving the "quadratic programming". Shorter code tends to be easier to write and maintain, assuming that readability is the same.

It is of course difficult to make a fair comparison of code lengths, since it is largely influenced by writing style. Still, there are a few things about tomSym that makes it inherently more suited for optimization than the symbolic toolbox.

For example, consider the following lines from the constraint definitions by DeLand:

```
% 50000 <= Stor(t) <= 100,000 AF
% Stor(t) = s0+sum(inFlow(1:t))-sum(spillFlow(1:t))-sum(turbineFlow(1:t))
c = stor0 + C2A*cumsum(inFlow); % Convert CFS to AF
b = [b; 100000-c; -50000+c];
s = -C2A*sparse(tril(ones(N)));
s = [s s];
A = [A; s; -s];
```

These five lines of code would be very hard to understand without the two lines of comments that accompany them. Even with the comments, it is not obvious how the constraint works.

With tomSym all of the above is replaced by a single line in the constraints array:

` 50000 <= storage <= 100000`

(This is an abuse of Matlab syntax. It should be written as two separate inequalities, but the ezsolve function accepts inequalities chained in this way, and it is a commonly used mathematical notation.)

This single line of code is clearly easier to read, write and debug than the seven lines above. Constraints can be written in this way because the ezsolve function automatically parses linear constraints and generates the necessary numeric matrices for the solver. In the same manner it also automatically generates automatic first and second derivatives of both objective and constraints, and passes these to the solver as matlab functions or numeric data depending on the problem type.

## Timing measurments

The symbolic toolbox needed a little over *5.1 hours* to compute the symbolic Hessian for 480 time steps (960 decision variables).
The quadprog solver then solved the numeric problem in about 10 seconds.

With TOMLAB, the symbolic Hessian was computed in 2.1 seconds, and the KNITRO numeric solver required 3.3 seconds, for a total
of 5.4 seconds. *An overall speedup of a factor 3 500.*

In to understand the difference between TOMLAB to the Mathworks' toolboxes, we ran a series of tests, using smaller problem sizes. Each test was run 5 times, and the mean of the middle 3 is listed in the table below. (N is the number of time steps. The number of decision variables is 2N.)

N | Symbolic TB | tomSym |
---|---|---|

10 | 0.698 s | 0.502 s |

50 | 12.86 s | 0.502 s |

100 | 67.49 s | 0.518 s |

150 | 237.0 s | 0.564 s |

The time for deriving the Hessian using the symbolic toolbox increased approximately as N^2.7 as the number of variables increased. With tomSym, on the other hand, the time remained approximately constant for problems this small.

N | quadprog | KNITRO |
---|---|---|

10 | 0.0162 s | 0.0165 s |

50 | 0.0452 s | 0.0318 s |

100 | 0.2015 s | 0.0829 s |

150 | 0.4756 s | 0.1847 s |

KNITRO runs at about the same speed as quadprog on the smallest problem sizes, and about 2.5 times faster for larger problems. (This is perhaps not a fair comparison, because quadprog has to do its own pre-scaling, while KNITRO gets a pre-scaled problem from tomSym.) In every case, the time for the numeric solver is much smaller than the time for symbolic processing.

## Comparing tomSym to the symbolic toolbox

The main reason why tomSym is so much faster than the symbolic toolbox is that it can handle vector and matrix symbols. This means that symbolic manipulations of matrix expressions are handled at approximately the same speed, regardless of the sizes of the symbolic matrices involved (until the matrices become so large that the time for numeric computations begins to matter.)

The symbolic toolbox, on the other hand, only handles scalar symbols. This means that both the number of symbolic expressions and their sizes increase as the size of symbolic matrices increase. The matlab code generated by the symbolic toolbox also becomes quite lengthy because of the use of multiple scalar operations instead of vector/matrix operations.

## Footnote

The above code passes the same Hessian to the solver as DeLand's code. It is possible to get a much sparser Hessian, by making the water level a decision variable instead of a symbolic expression. With tomSym it is very easy to to this kind of modification. All we need to do is to replace the expression for "storage" by a symbol definition:

` storage = tom('storage',N,1);`

and then add one constraint

` constraints{end+1} = ( storage == [stor0; storage(1:N-1)] + inFlow - outFlow );`

This modification increases the number of decision variables from 2N to 3N, but still makes both the symbolic processsing and the numeric solver run much faster.

For N=3200 (corresponding to ca 4.4 months of data), a solution was obtained 6.4 seconds. (5.4 seconds symbolic processing and 1 second for the KNITRO numeric solver.)