Constrained dynamic optimization with fmincon

Hello,

I am trying to use fmincon to solve a dynamic optimization problem:

where u are my control variables, x are my state variables and T are terminal constraints that depend on the state variables. I know that in theory in a dynamic optimization problem there are adjoint variables that vary with time and that are defined with regards to the constraints. The optimality conditions depend on those adjoint variables (optimality conditions based on Pontryagin's principle).

In my case, I need T to be equal to 0 (within a specific tolerance). Not knowing how to solve a dynamic optimization problem with Matlab, I have taken these contraints into consideration by creating a new optimization function Jnew= w1*J+ w2 * T and I try to minimize Jnew.

For example I set w2 by keeping the tolerance I want for T in mind, e.g. if I want T to be around 1e-3 and w1*J has an order of magnitude of 1 at optimum (based on calculations run previously without the contraint), I set w2 to 1e3 to have the same order of magnitude for my two terms at the "optimum", so that fmincon doesn't just consider one term over the other. However I feel this is a somewhat arbitrary way to set the weights. And the problem is that they remain constant in all iterations.

I just don't know if this is a legitimate way to apply fmincon, if the way I have formulated my problem allows me to obtain a minimum that satisfies the optimality conditions. I know fmincon is not a dynamic solver and it doesn't consider adjoint variables, especially since I directly define the constraints myself as a sum.

Can anyone with some knowledge on the subject or has been confronted with a similar problem help me out please?

I would really appreciate it.

Answer 1

[Torsten]

J is minimized in the objective function.

The constraint T(x(t_final)) = 0 is set in the function "nonlcon" of fmincon.

u (and if there are constraints on x, also x) is/are the vector(s) of unknowns to be determined by fmincon.

Constraints on u and/or x can be set in lb, ub, A, b, Aeq,beq and nonlcon (depending on their nature).

So don't mix objective function and constraint by building O = c1*J + c2*T as objective - J is the objective and T=0 is a constraint.

[AM]

Hi,

x are state variables determined by integration of the state equations dx = f(x,u,t). I use ode15s to solve the state equations. The state variables are not to determined by fmincon.

As I mentioned, it is a dynamic optimization problem.

[Torsten]

If you want to put constraints on the state variables, you will have to treat them as variables in "fmincon". This is usually done by imposing constraints on x at a grid of intermediate integration times: x(t1),x(t2),...,x(tf). "fmincon" will then try to adjust the control variables u such that the constraints on x(t1), x(t2),...,x(tf) are satisfied.

If not, you don't have to put x (at discrete times ti) in the list of "fmincon" variables.

From your problem description it seems that the only constraint is T(x(tf)) = 0 and that the x-variable(s) are unconstrained. But this is rather unusual.

[AM]

Page 35 in this document presents the type of problem I have: Dynamic optimization with constraints. Maybe it's more clearly presented here.

If I add the state variables as variables in fmincon won't fmincon try to make them vary on its own as its does its integrations? I just don't see how I can add them as optimization variables but let my code know they are not to be calculated by the optimization algorithm but rather by integration of the state equations. Do you have an example showing how to do what you suggest?

Thank you for your help by the way.

[Torsten]

Given control parameters u(ti), in nonlcon, you integrate the system xdot = f(x,u,t), x(t0) = x0 independent of the variables x_given_from_fmincon(ti).

After integration, you supply the equality constraints

x(ti) - x_given_from_fmincon(ti) = 0

and maybe inequality constraints

something1(ti) <= x_given_from_fmincon(ti) <= something2(ti)

These two settings force fmincon to adjust the control variables u(ti) (if possible) such that the integration gives

something1(ti) <= x(ti) <= something2(ti)

But - as said - you don't need to define the x(ti) as unknowns for "fmincon" if you don't want to impose constraints on them.

[AM]

I've applied it what you suggested to a simple example and I've managed to make it work, thank you !

I just have one more question. My real problem is very time consuming because of the integration of the state variables and I was wondering if the differential system is solved twice with this method, as it is solved within the function that calculates the objective function (objfunction in the example below) and then it is solved within the nonlcon function.

Essentially, I have two functions that both integrate the differential system, is this the correct way to do it?

function [c,ceq] = nonlcon(controlvector,extraParams)

[t,y] = ode15s(@odefun,tspan,y0); %tspan and y0 are included in the extraParams vector

c = []

ceq = y(end)-valuecons; % I set y at tf to be equal to valuecons

end

function optvariable = objfunction(controlvector,extraParams)

[t,y] = ode15s(@odefun,tspan,y0); %tspan and y0 are included in the extraParams vector

optvariable=f(y) % my optimization variable is calculated by using the values of y

end

[Torsten]

Yes, the differential system must be solved twice - once in "nonlcon" and once in "objfunction".

And if you now want to ask whether you can save one call to ode15s, my answer is: contact MATLAB support because it can be true that both functions are called in a certain order and frequency with the same input arguments, but I don't know the order and the frequency.

What would be a good programming style is to write a function "state_space_integration" that performs the integration and is called both from "nonlcon" and "objfunction" to get y. But this will not save computing time, of course.

[Walter Roberson]

Use memorization to avoid extra computation.

[Torsten]

@Walter Roberson

That means saving p and y of the last calls to the "state_space_integration" function and comparing to p of the actual call ?

[AM]

I found this when I looked for a way to solve my problem : Objective and Nonlinear Constraints in the Same Function. I've implemented it and it does the job perfectly.