TomSym compared to the Symbolic Toolbox: Difference between revisions
No edit summary |
|||
(41 intermediate revisions by the same user not shown) | |||
Line 3: | Line 3: | ||
We look at the optimization problem of maximizing the revenue of a hydroelectric dam. | 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. [http://www.mathworks.com/company/newsletters/articles/solving-large-scale-optimization-problems-with-matlab-a-hydroelectric-flow-example.html] | This problem was published by By Seth DeLand of the MathWorks, as an example of how to combine their optimization and symbolic toolboxes. [http://www.mathworks.com/company/newsletters/articles/solving-large-scale-optimization-problems-with-matlab-a-hydroelectric-flow-example.html] | ||
We solve the same probem using tomSym and the TOMLAB solvers, noting a significant speed difference | 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 == | == Problem description == | ||
Line 21: | Line 21: | ||
* The final reservoir level must equal the initial level (90000 AF) | * The final reservoir level must equal the initial level (90000 AF) | ||
The | The code and data used by DeLand can be downloaded from Matlab Central.[http://www.mathworks.com/matlabcentral/fileexchange/35856-optimization-in-matlab-an-introduction-to-quadratic-programming] | ||
== tomSym code == | == tomSym code == | ||
Line 55: | Line 55: | ||
electricity = turbineFlow.*(0.5*k1*(storage + [stor0; storage(1:end-1)]) + k2)./MW2kW; | electricity = turbineFlow.*(0.5*k1*(storage + [stor0; storage(1:end-1)]) + k2)./MW2kW; | ||
revenue = price*electricity; | revenue = price'*electricity; | ||
objective = -revenue; % Minimization objective | objective = -revenue; % Minimization objective | ||
Line 71: | Line 71: | ||
options = struct; | options = struct; | ||
options.solver = 'snopt'; | options.solver = 'knitro'; % also try: 'snopt' | ||
options.Prob.KNITRO.options.ALG=3; | |||
options.scale = 'auto' | options.scale = 'auto' | ||
Line 81: | Line 82: | ||
[[File:Hydroelectric plots a.png|thumb|Plots of the result, when the approximate indata is used.]] | [[File:Hydroelectric plots a.png|thumb|Plots of the result, when the approximate indata is used.]] | ||
Although the Hessian is | Although the Hessian is constant, this is <em>not</em> a quadratic programming. | ||
The Hessian is not positive semidefinite, and the objective is therefore not convex. There may exist | 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 | local optima that are not as good as the global one. For this reason, the QP solvers in the | ||
Line 87: | Line 88: | ||
The solution for the approximate input data is very similar to the one presented by DeLand. | 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 | 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 | |||
<source lang="matlab"> | |||
(subs(Rtot,X,x3)-fval3)./fval3 | |||
</source> | |||
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 == | == Code comparison == | ||
Line 94: | Line 103: | ||
One thing that is immediately noticeable when comparing the above code to DeLand's is that it is much shorter. | 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 | 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". | 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. | Shorter code tends to be easier to write and maintain, assuming that readability is the same. | ||
Line 112: | Line 121: | ||
</source> | </source> | ||
These five lines of code would be very hard to understand without the two lines of comments that accompany them. | 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 | 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: | With tomSym all of the above is replaced by a single line in the constraints array: | ||
Line 118: | Line 127: | ||
50000 <= storage <= 100000 | 50000 <= storage <= 100000 | ||
</source> | </source> | ||
(This is an abuse of Matlab syntax. It should be written as two separate | (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 == | == Timing measurments == | ||
In | The symbolic toolbox needed a little over ''5.1 hours'' to compute the symbolic Hessian for 480 time steps (960 decision variables). | ||
Each test was run 5 times, and the mean of the middle 3 | 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.) | |||
{| class="wikitable" | {| class="wikitable" | ||
|+ Time for symbolic processing | |+ Time for symbolic processing | ||
|- | |- | ||
! | ! N || Symbolic TB !! tomSym | ||
|- | |- | ||
| 10 || 0.698 s || 0.502 s | | 10 || 0.698 s || 0.502 s | ||
Line 141: | Line 158: | ||
The time for deriving the Hessian using the symbolic toolbox increased approximately as N^2.7 as the number of variables increased. | 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. | With tomSym, on the other hand, the time remained approximately constant for problems this small. | ||
{| class="wikitable" | {| class="wikitable" | ||
|+ Time for numeric solver | |+ Time for numeric solver | ||
|- | |- | ||
! | ! N || quadprog !! KNITRO | ||
|- | |- | ||
| 10 || 0.0162 s || 0.0165 s | | 10 || 0.0162 s || 0.0165 s | ||
Line 157: | Line 174: | ||
|} | |} | ||
KNITRO | 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 == | == 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. | 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 matrix expressions are handled at approximately the same speed, regardless of the sizes | 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.) | ||
of the symbolic matrices involved. | |||
The symbolic toolbox, on the other hand, only handles scalar symbols. This means that both the number of | The symbolic toolbox, on the other hand, only handles scalar symbols. This means that both the number of | ||
expressions and their sizes increase as the size of symbolic matrices increase. | 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 | 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. | scalar operations instead of vector/matrix operations. | ||
Line 174: | Line 190: | ||
The above code passes the same Hessian to the solver as DeLand's code. | 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. | It is possible to get a much sparser Hessian, by making the water level a decision variable instead of a symbolic expression. | ||
All we need to do is to replace the | 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: | |||
<source lang="matlab"> | <source lang="matlab"> | ||
storage = tom('storage',N,1); | storage = tom('storage',N,1); | ||
Line 182: | Line 199: | ||
constraints{end+1} = ( storage == [stor0; storage(1:N-1)] + inFlow - outFlow ); | constraints{end+1} = ( storage == [stor0; storage(1:N-1)] + inFlow - outFlow ); | ||
</source> | </source> | ||
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.) | |||
[[Category:tomSym]] |
Latest revision as of 08:34, 21 January 2015
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.)