Optimization: Reuse calculations from objective function in nonlinear constraints function

Hello,

Is it possible to reuse calculations from one function in the other to save computational effort?

Example:

function SSE = objfunc(inputs)
% some calculations that output model results + some other things
SSE = sum(residuals.^2)
end

function [c, ceq] = nonlinconstraint(inputs)
% same calculations as the objective function
building c and ceq vectors
end

Now, I know I could make a third function that wraps the calculations I mention and call that function both in objfunc and nonlinconstraint but that would imply in doing (the same) expensive calculations twice. Is there a way to reuse the calculations from one function into another or calculate it only once and pass the results to them? Other some other workflow that saves time?

If I use parallel computing toolbox I doubt it would reuse these common calculations, although it would probably make the code faster.

I am talking in the context of optimization, like multistart, globalsearch, fmincon, etc. that take up those functions to try to find global/local optimum.

I could potentially write a single function that outputs SSE,c,ceq but then I do not know how I would pass it to the aforementioned algorithms.

Answer 1

[Torsten]

https://de.mathworks.com/help/optim/ug/objective-and-nonlinear-constraints-in-the-same-function.html

[Gustavo Lunardon]

Thank you Torsten! This is much in line with what I was looking for, glad to see there is a documented example on that. I was not aware. I will try to make a mix of this example and Bruno's suggestion to make it work for my case.

[Gustavo Lunardon]

Does this also work with linear constraints or the solvers call the linear constraints more frequently than objective/nonlinear constraints?

[Torsten]

Linear constraints are set in A,b,Aeq,beq,lb and ub.

Since the arrays/matrices are constant, they don't have to be updated.

[Gustavo Lunardon]

OK. I see I can construct things that way. Thanks

[Gustavo Lunardon]

Actually, I have equality constraints in the form Aeq*x = beq(x). I was using the nonlinear constraints to pass that to the solver. I see that fmincon and the other algorithms do not accept function handles for beq so guess there is no other way around it.

[Torsten]

Yes, Aeq*x = beq(x) is a nonlinear constraint to be defined in "nonlinconstraint".

But why do you think you could save time here ?

You can't reuse calculations already performed in "objfunc" in this case, do you ?

[Gustavo Lunardon]

I do, the constraints are in that form because I am trying multiple shooting. They are continuity constraints and I need to integrate from the starting point of the node. These calculations yield the model value at the point (used for the residual in the objective function) and the values for the states at each node, which should be equal to the final value of the integration of the previous node (hence the nonlinear constraint). It does save time if I do not need to integrate everything twice.

Well, I could pass the initial and final state constraints as linear constraints, even if I could not do the same with the continuity ones.

[Torsten]

And you think multiple shooting is more stable than just simple integration using one of the ODE integrators ? Or why did you choose such a complicated approach ?

[Gustavo Lunardon]

I can use one of the standard ODE integrators between the nodes, no need to code Runge-Kutta myself, although I have that too. I do not think multiple shooting is more stable, I am struggling to make it converge at the moment. I am doing this because I am trying to solve optimal control problems and the control variables are discretized in time. I was trying to apply multiple shooting for regression to familiarize myself with it.

So control problems of the form, just as an example:

dy(1) = y(1);

dy(2) = u;

u could be a sequence of piecewise constant accelerations and y(1) and y(2) position and velocity, respectively. Then minimize some function related to the fuel used and so on.

I could also try to code orthogonal collocation but that looks more complicated than multiple shooting.

[Torsten]

But a usual ODE integrator also accepts discretized inputs. You just have to define how they should be interpolated between the discretization points. But there are several options: linear, quadratic, spline ...

But you must know the way to go with your problem.

[Gustavo Lunardon]

Do you mean define a vector u = [u1,u2,u3,...,uN], where N = number of discretrizations for the control variable, pass it as an additional parameter to the ode function and select (via interp1, pchip, etc) which element of the vector to be used in the integration for that given time step (internal time step as selected automatically by the ode solver)? Or something else?

[Torsten]

You get the control vector

u = [u1,u2,u3,...,uN]

passed from "fmincon", I guess, because it is usually made up of variables to be adapted/optimized. With this control vector, one usually forms a function handle

fun_u = @(t)interp1(tspan,u,t,'spline') (or something else instead of spline)

passes this function handle to the ode-integrator (more precisely: to the function where you define the ode system) and evaluates it at the places where the control variable needs to be inserted as

u = fun_u(t).

Note that the (length of u) might be the (length of tspan - 1) ; this has somhow to be handled in the definition of fun_u.

[Gustavo Lunardon]

Thanks for the explanation! Passing a handle looks better than passing a vector as I was thinking. Looks more neat.

Answer 2

[Bruno Luong]

You can use this design pattern

function preciousvalue = expensivereuse(inputs)
persistent lastresult
if ~isempty(lastresult) && isequal(input, lastresult.inputs)
    preciousvalue = lastresult.preciousvalue;
else
    preciousvalue = computepreciousvalue(inputs); % this is the function that computes the expensive part
    lastresult = struct('inputs', inputs, 'preciousvalue', preciousvalue);
end
end

function SSE = objfunc(inputs)
preciousvalue = expensivereuse(inputs);
% some specific calculations that output model results + some other things
% from preciousvalue
SSE = sum(residuals.^2)
end

function [c, ceq] = nonlinconstraint(inputs)
preciousvalue = expensivereuse(inputs);
% some specific calculations 
from preciousvalue
building c and ceq vectors
end

[Gustavo Lunardon]

That is an interesting approach. It uses some lines I was not familiar with, like the persistent and from. I will keep that in mind. Thanks for the suggestion! It opens some paths to do other things I was looking for as well.