# Why is optimization involving ODEs so much slower than with a purely algebraic objective function?

1 view (last 30 days)
rotton on 25 Feb 2019
Commented: dpb on 25 Feb 2019
Consider this optimization script, which is a very much boiled down version of something I am currently working on:
% Parameters
substr_0 = 6.21;
substr_in = 3.171;
k_1 = 0.5;
f_1 = 43;
phi_1 = 0.5;
initialVolume = 0;
daysTotal = 7;
hoursTotal = daysTotal*24;
%% Consumer schedule per hour
onOff = [0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 ...
1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 1 1 1 1 1 ...
1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 ...
1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ...
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0]';
% Initial value: Feeding only during "Consumer ON" time
substr_in_time = substr_in/12.*onOff;
%% Setup optimization problem
parameters = [substr_0, k_1, f_1, phi_1, initialVolume];
% Objective function: Standard deviation of storage volume
objFun = @(feedingVector) substr_Model(feedingVector, onOff, parameters);
% Bounds (lb <= x <= ub), i.e., decision variables can only range between 0 and 6.3
lowerBound = zeros(1, hoursTotal);
upperBound = 2*substr_in.*ones(1, hoursTotal);
% Possible solvers: FilterSD, Ipopt, fmincon (local)
optimOpts = optimoptions('fmincon', 'Display', 'iter', 'FunctionTolerance', 1e-2,...
'MaxIterations', 1e5, 'MaxFunctionEvaluations', 1e5*daysTotal);
%% Solve optimization problem
tic; [feeding_Opt, fval] = fmincon(objFun, substr_in_time, [], [], [], [],...
lowerBound, upperBound, [], optimOpts); toc
%% Calculate model with optimal values
substr_Model(feeding_Opt, onOff, parameters)
%% Plot results
figure();
plot(1:hoursTotal, interm1, '-o', 'DisplayName', 'Interm production');
hold on; grid on
stairs(0:hoursTotal-1, 5*onOff, '-', 'DisplayName', 'Consumption schedule');
plot(1:hoursTotal, storageVolume/10, '-^', 'MarkerSize', 4,...
'DisplayName', 'Storage volume [1e-1*m^3]');
stairs(0:hoursTotal-1, feeding_Opt, 'DisplayName', 'Feeding schedule');
xtickVector = 0:12:hoursTotal;
xticks(xtickVector);
xticklabels(xtickVector/24);
xlabel('Time (days)');
legend('Location', 'Best');
set(gca, 'FontSize', 18);
hold off
function storageStd = substr_Model(substr_feeding, onOff, parameters)
hoursTotal = length(substr_feeding);
substr_0 = parameters(1);
k_1 = parameters(2);
f_1 = parameters(3);
phi_1 = parameters(4);
initialVolume = parameters(5);
% Simple for-loop with feeding as vector per hour
substr = zeros(hoursTotal, 1);
substr(1) = (substr_0 + substr_feeding(1))*exp(-k_1/24);
for hour = 2:hoursTotal
substr(hour) = (substr(hour-1) + substr_feeding(hour))*exp(-k_1/24);
end
% Calculate resulting interm1 production and storage volume
interm1 = phi_1*f_1*substr/24;
interm1ConsumptionRate = 11.24509;
storageVolume = initialVolume + cumsum(interm1 - onOff*interm1ConsumptionRate);
assignin('base', 'substr', substr);
assignin('base', 'interm1', interm1);
assignin('base', 'storageVolume', storageVolume);
storageStd = std(storageVolume);
end
However, If I use an ODE as part of the objective function, optimization takes much longer!
function storageStd = substr_Model(substr_feeding, onOff, parameters)
hoursTotal = length(substr_feeding);
substr_0 = parameters(1);
k_1 = parameters(2);
f_1 = parameters(3);
phi_1 = parameters(4);
initialVolume = parameters(5);
% Simple for-loop with feeding as vector per hour
substr = zeros(hoursTotal, 1);
initialValues4Model_0 = [substr_0 + substr_feeding(1)];
% Solve ODE
[~, x_sol] = ode15s(@(t, x) substr_Model_diff(t, x, parameters), 0:0.5:1,...
initialValues4Model_0);
substr(1) = x_sol(end);
for hour = 2:hoursTotal
initialValues4Model = [substr(hour-1) + substr_feeding(hour)];
[~, x_sol] = ode15s(@(t, x) substr_Model_diff(t, x, parameters), 0:0.5:1,...
initialValues4Model);
substr(hour) = x_sol(end);
end
% Calculate resulting interm1 production and storage volume
interm1 = phi_1*f_1*substr/24;
interm1ConsumptionRate = 11.24509;
storageVolume = initialVolume + cumsum(interm1 - onOff*interm1ConsumptionRate);
assignin('base', 'substr', substr);
assignin('base', 'interm1', interm1);
assignin('base', 'storageVolume', storageVolume);
storageStd = std(storageVolume);
end
With the ODE in a separate file:
function dx = substr_Model_diff(t, x, parameters)
k_1 = parameters(2);
% RHS of ODE
dx = -k_1/24*x;
end
(Differential objective code could probably be improved a bit, but this is just for demonstration purposes)
So my question is: Why is optimization involving an ODE so much slower than with a purely algebraic objective function? I understand fmincon uses the obejctive's Jacobian. Is it not able to determine this automatically in the second case? Would passing an analytic Jacobian (created from Symbolic Math Toolbox) help? But then, how to derive it, since the objective contains a for loop?
dpb on 25 Feb 2019
How could it avoid being "so much slower" when the ODE has to be solved to get the functional value to optimize over?

R2018a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!