How to constrain a vector when using fmincon?

I wrote a fmincon to calculate the optimal inputs:

[ In_opt ] = fmincon( @( Inputs_vector )obj( other parameters_1, Inputs_vector), [Inputs_vector_0], [], [], [], [], [Inputs_vector_L],[Inputs_vector_U],... @(  )ctr(  ), optNLP);

There is an integration process when calculating function "obj", for example, t = 0:0.02:2, I used

    for

    ode45

    end

during this process. For every step of ode45, for example, from 0.02s to 0.04s, there is an intermediate variable K (depending on current system states and input) which I want to constrain between [-3, 3]. Here my total simulation time is 2s, thus for the simulation from 0 to 2s, there are about 101 K values.

I want to make sure that when I use fmincon, for every time interval: [0, 0.02], [0.02, 0.04], [0.04, 0.06]......[1.96, 1.98], [1.98, 2], the K are constrained within [-3, 3]. But how should I write this constrain condition in @( )ctr( )?

• "ctr" is the constrain function for my fmincon.

Thank you very much, I'm waiting for your reply.

Answer 1

[Alan Weiss]

I think that you need to do two things:

1. Decide how many points on the ODE solution that you want to constrain
2. Use the documented method to run your constraints and objective function in the same call

Suppose that you decide, as you suggested, that you want to constrain 100 time steps. The, when you run the ODE solver, extract the solution at the 100 points. And create 200 nonlinear constraints. One nonlinear constraint is the vector of lower bounds at each point:

 -3 - K(t)

Here I am assuming that K is a 100-long vector of the K values. As long as each K(t) is greater than -3 then -3-K(t) < 0, which is what the nonlinear constraint function wants. But if some K(t) < -3 then -3-K(t) > 0, which means a nonlinear inequality constraint violation.

For the upper bounds, use

 K(t) - 3

As long as K(t) < 3 then K(t) - 3 < 0.

And to save time, be sure to use the method for calling constraints and objective in the same function call.

Alan Weiss

MATLAB mathematical toolbox documentation

[Alan Weiss]

I should also add that you might need to take larger-than-default gradient estimation steps, as explained in Optimizing a Simulation or ODE.

Answer 2

[Brendan Hamm]

Please format your question so it is easier to read.

Is K one of the design variables for the fmincon call? If so, then you just need to set the lower bounds and upper bounds in fmincon. The function nonlcon passed to the fmincon routine is just for non-linear constraints.

fun = @( Inputs_vector )obj( other parameters_1, Inputs_vector);

x0 = Inputs_vector_0;

lb = -3*ones(length(K));
ub = 3*ones(length(K));

x = fmincon(fun,x0,[],[],[],[],lb,ub,...) % Complete with other inputs

If K is not a design variable, then this question needs to be reworded/re-titled.

[Tao]

K is not designed, but a variable always depending on current step system states, and K is changing when solving the ode from 0s to 2s. For example:

K(i) = (X(1) + X(2)*a)/X(1);

X(1) is state 1 (scalar), X(2) is state 2 (scalar). i is current ode step (from 0s to 2s, there are 100 steps to integrate the system, so X(1) and X(2) are always changing with step too)

I want to constrain K in every step of ode through optimize the input which will influence X(1) and X(2).

[Brendan Hamm]

This is still not very clear. Can you please post your code?

[Tao]

First, thank you for your patience. The original code is long and complicated, I will write a simplified one here, but it may be not logically accurate, just to show what I mean:

The main code:

Time = 0:0.02:2;

[ In_opt ] = fmincon( @( Inputs_vector )obj( other parameters_1, Inputs_vector), [Inputs_vector_0], [], [], [], [], [Inputs_vector_L],[Inputs_vector_U],... @( )ctr( ), optNLP);

Objective function:

function [ J ] = obj( other parameters_1, Inputs_vector )

[ f ] = fun( other parameters_1, Inputs_vector );

J = -f(1);

end

function "fun.m": (which is used to do the integration from 0s to 2s)

function [ f ] = fun( other parameters_1, Inputs_vector )

for i = 1:1:(length(Time)-1)

X(2) = g(X(1), Inputs_vector); % g is a complicated nonlinear function

K(i) = (X(1) + X(2)*a)/X(1);

Tspan = [Time(i) Time(i+1)];

[~, X] = ode45(@(t,X) state( parameters_decided_by_states_and_inputs ), Tspan, IC);

X(1) = X(end, 1);

X(2) = X(end, 2);

X(3) = X(end, 3);

end

f(1) = X(1);

end

The "state" is some differential equations (with time).

Answer 3

[Matt J]

If you want K bounded at a discrete set of times, just generate the stream of K at those times inside ctr and set

    cineq=[K(:)-3; -3-K(:)]

If K has to be bounded over a continuum of times, you may need to use fseminf().

[Tao]

Could you elaborate more?

What I want is, for the simulation process time intervals:

[0, 0.02], [0.02, 0.04], [0.04, 0.06]......[1.96, 1.98], [1.98, 2], all k in these intervals are within [-3, 3]. Since k is influence by system states and system control input, that's why I use fmincon to get the optimal system control input.

[Matt J]

Just as you've generated the K(i) as a function of x inside your objective function, you can generate them as well inside your nonlinear constraint function. Once you do that, you can build the inequality constraint vector cineq appropriately based on the vector K. For example,

   function [cineq,ceq]=ctr(x)

     ...
     ceq=[];

     cineq=K-3;

   end

is equivalent to saying that all K(i)<=3. Although, as Alan also mentioned, you can avoid computing the K vector twice by using techniques described here.