improving symbolic generation of objective function for optimization

Question

I am currently using matlab to solve an optimization problem. I am generating the objective function using the symbolic toolbox. I planned use the symbolic toolbox to calculate the gradient and hessian and so speed up the optimization. The problem I have however is that when N is large (in my case over 8000) generating the function handles takes hours.

So my questions are these:

• is there a way to speed up function handle generation?
• is there an alternate way of generating an objective function when the length of N varies and be able to find the gradient and hessian?
• if not perhaps I should be using another package recommended for such a task? if so any recommendations?

here is a code snippet to show what I mean

X1 = sym('x1',[N,1]);
X2 = sym('intx',[2*N,1]);
P = sym('P',[N,1]);
FC = sym('FC',[N,1]);
CC = sym('CC',[N,1]);
SC = sym('SC',[N,1]);
AC = sym('AC',[N,1]);

efficiency6 = 0.1199 * (X1(1:N)/max_el_capacity_6).^3 - 0.3568 * (X1(1:N)/max_el_capacity_6).^2 + 0.4031* (X1(1:N)/max_el_capacity_6) + 0.2286;

income6 = X1(1:N).*(P(1:N)-CC(1:N)-AC(1:N)-FC(1:N)./efficiency6(1:N));

revenue6 = X2(1:N).*(income6) - SC(1:N).*X2(N+1:2*N);

totrevenue6 = -sum(revenue6);

totRevenue6 = subs(totrevenue6,[P;CC;FC;AC;SC],[ep';cc';fc';ac6';sc6']);

matlabFunction(totIncome6,'vars',{X1},'file','objectiveFcn2014_1');

many thanks, Jesse

Answer 1

Using symbolics to derive the hessian of a problem with 8000 variables sounds like a bad idea, unless it is extremely sparse. For problems of that size, if you absolutely want to work with the exact hessian, you would go for a tool which performs automatic differentiation (which is something different from difference approximations and symbolic computation of the whole hessian). I only know of ADImat along that direction.

If you simply want to have a convenient modelling language inside MATLAB with connections to nonlinear solvers (I let other guide you outside MATLAB), you could go for YALMIP (disclaimer, developed by me). Derivatives will be supplied automatically to the solver, but the solver will have to rely on quasi-newton methods for the hessian

X1 = sdpvar(N,1);
X2 = sdpvar(2*N,1);
P = sdpvar(N,1);
...

efficiency6 = 0.1199 * (X1/max_el_capacity_6).^3-0.3568 * (X1/max_el_capacity_6).^2 +     0.4031* (X1/max_el_capacity_6) + 0.2286;
...
% maximize revenue6 s.t. P >=0
% fmincon most likely used if that is what you have installed
Constraints = [P>=0];
Objective = -sum(revenue6);
optimize(Constraints,Objective)

If you want to pursue this, you can contact me in private for a detailed discussion which might be uninteresting on this thread.