how to solve a multi-objective nonlinear optimization problem with constraints ?

Hello , I have a nonlinear eqt that describes the evolution of my system : dxdt=Ax(t)+Bu(t),with

A=[sqrt(2) 1;2 sqrt(3)]; B=[1 4;3 1]; and u=K1*x+b1. K1 is 2x2 matrix and b1 is 2x1.

I want to find K1 and b1 such that I maximize xdot11 at a point v1 under these constraints :

n1'*B*(K1*v1+b1)+n1'*A*v1<=0;

n4'*B*(K1*v1+b1)+n4'*A*v1<=0;

n1'*xdot11<= 0 ;

K1*v1+b1 <= 2;

K1*v1+b1 >=- 2;

n4'*xdot12<= 0 ;

n1=[0;-1]; n4=[-1;0]

My code is :

clc;
close all;
v1=[80;60]; 
n1=[0;-1]; n4=[-1;0]; % normal vectors
A=[sqrt(2) 1;2 sqrt(3)];
B=[1 4;3 1];

prob1 = optimproblem;
K1 = optimvar("K1",2,2); % create the optim. variable K1
b1 = optimvar("b1",2,1);
xdot11=((A(1,1)+B(1,1)*K1(1,1)+B(1,2)*K1(2,1))*v1(1))+((A(1,2)+B(1,1)*K1(1,2)+B(1,2)*K1(2,2))*v1(2))+B(1, ...
  1)*b1(1)+B(1,2)*b1(2);
xdot21=((A(2,1)+B(2,1)*K1(1,1)+B(2,2)*K1(2,1))*v1(1))+((A(2,2)+B(2,1)*K1(1,2)+B(2,2)*K1(2,2))*v1(2))+B(2, ...
    1)*b1(1)+B(2,2)*b1(2);

Obj =-xdot11;
prob1.Objective=Obj

prob1 = 
  OptimizationProblem with properties:

       Description: ''
    ObjectiveSense: 'minimize'
         Variables: [1×1 struct] containing 2 OptimizationVariables
         Objective: [1×1 OptimizationExpression]
       Constraints: [0×0 struct] containing 0 OptimizationConstraints

  See problem formulation with show.

constr1=n1'*B*(K1*v1+b1)+n1'*A*v1<=0;
constr2=n4'*B*(K1*v1+b1)+n4'*A*v1<=0;
constr3=n1'*xdot11<= 0 ;
constr4=K1*v1+b1 <= 2;
constr5=K1*v1+b1 >=- 2;
constr6=n4'*xdot21<= 0 ;
prob1.Constraints.Constr1=constr1;
prob1.Constraints.Constr2=constr2;
prob1.Constraints.Constr3=constr3;
prob1.Constraints.Constr4=constr4;
prob1.Constraints.Constr5=constr5;
prob1.Constraints.Constr6=constr6;

show(prob1);

  OptimizationProblem : 

	Solve for:
       K1, b1

	minimize :
       -80*K1(1, 1) - 320*K1(2, 1) - 60*K1(1, 2) - 240*K1(2, 2) - b1(1) - 4*b1(2) - 173.1371


	subject to Constr1:
       -240*K1(1, 1) - 80*K1(2, 1) - 180*K1(1, 2) - 60*K1(2, 2) - 3*b1(1) - b1(2) <= 263.923

	subject to Constr2:
       -80*K1(1, 1) - 320*K1(2, 1) - 60*K1(1, 2) - 240*K1(2, 2) - b1(1) - 4*b1(2) <= 173.1371

	subject to Constr3:
       -80*K1(1, 1) - 320*K1(2, 1) - 60*K1(1, 2) - 240*K1(2, 2) - b1(1) - 4*b1(2) <= 173.1371

	subject to Constr4:
       80*K1(1, 1) + 60*K1(1, 2) + b1(1) <= 2
       80*K1(2, 1) + 60*K1(2, 2) + b1(2) <= 2

	subject to Constr5:
       80*K1(1, 1) + 60*K1(1, 2) + b1(1) >= -2
       80*K1(2, 1) + 60*K1(2, 2) + b1(2) >= -2

	subject to Constr6:
       -240*K1(1, 1) - 80*K1(2, 1) - 180*K1(1, 2) - 60*K1(2, 2) - 3*b1(1) - b1(2) <= 263.923
     

x0.K1=zeros(2); 
x0.b1=[1;1];
 sol =fsolve(prob1,x0); 

Error using fsolve
FSOLVE requires the following inputs to be of data type double: 'X0'.

when I run my code , I get this error :

Unable to perform assignment because dot indexing is not supported for variables of this type.

I tried to solve it with fmincon but it requires to put all of the control variables into one vector x and I didn't know how to define the x0 in this case..

I will be very thankfull for your help.

[Torsten]

You must set up the problem with the solver-based approach for fmincon.

How should MATLAB know from your settings above that it should solve a differential equation somewhere ?

[tty]

Thank you for your reply .I tried to solve it with fmincon but i don't know how shoud I define the intial value x0. Sorry if it is a dummy qst but this is my first time to solve a prob with fmincon.

function fonc(u1)
K1=[u1(1),u1(2);u1(3),u1(4)];
b1=[u1(5);u1(6)];
v1=[80;60];
A=[sqrt(2) 1;2 sqrt(3)];
B=[12.5 4;5 10];
xdot11=-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6));

end
function[c1,c2,c3,c4,c5,c6,ceq]=carreconst(u1)
v1=[80;60]; 
A=[sqrt(2) 1;2 sqrt(3)];
B=[1 4;3 1];
K1=[u1(1),u1(2);u1(3),u1(4)];
b1=[u1(5);u1(6)]';
xdot11=-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6));
xdot21=((A(2,1)+B(2,1)*u1(1)+B(2,2)*u1(3))*v1(1))+((A(2,2)+B(2,1)*u1(2)+B(2,2)*u1(4))*v1(2))+B(2, ...
    1)*u1(5)+B(2,2)*u1(6);
n1=[0;-1]; n4=[-1;0];

c1=n1'*B*(K1*v1+b1)+n1'*A*v1;
c2=n4'*B*(K1*v1+b1)+n4'*A*v1;
c3=n1'*xdot11;
c4=n4'*xdot21 ;
c5=K1*v1+b1-5;
c6=-K1*v1-b1+5;
ceq=[];

end

clc;
close all;
v1=[80;60];
n1=[0;-1]; n4=[-1;0]; % normal vectors
A=[sqrt(2) 1;2 sqrt(3)];
B=[1 4;3 1];
nonlcon=@carreconst;
fun=@(u1)(-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6)));
%x0=...? 
sol=fmincon(fun,x0,[],[],[],[],[],[],nonlcon)

[Torsten]

The expressions

n1'*B*(K1*v1+b1)+n1'*A*v1
n4'*B*(K1*v1+b1)+n4'*A*v1
n1'*xdot11
n4'*xdot21 
K1*v1+b1-5
-K1*v1-b1+5

you want to constrain are 1x2 vectors and 2x2 matrices. Is this really what you want ?

I still don't see that you solve your differential equation somewhere.

clc;
close all;
v1=[80;60];
n1=[0;-1]; n4=[-1;0]; % normal vectors
A=[sqrt(2) 1;2 sqrt(3)];
B=[1 4;3 1];
nonlcon=@carreconst;
fun=@(u1)(-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6)));
x0=ones(6,1); 
sol=fmincon(fun,x0,[],[],[],[],[],[],nonlcon)

function fonc(u1)
K1=[u1(1),u1(2);u1(3),u1(4)];
b1=[u1(5);u1(6)];
v1=[80;60];
A=[sqrt(2) 1;2 sqrt(3)];
B=[12.5 4;5 10];
xdot11=-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6));

end
function[c,ceq]=carreconst(u1)
v1=[80;60]; 
A=[sqrt(2) 1;2 sqrt(3)];
B=[1 4;3 1];
K1=[u1(1),u1(2);u1(3),u1(4)];
b1=[u1(5);u1(6)]';
xdot11=-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6));
xdot21=((A(2,1)+B(2,1)*u1(1)+B(2,2)*u1(3))*v1(1))+((A(2,2)+B(2,1)*u1(2)+B(2,2)*u1(4))*v1(2))+B(2, ...
    1)*u1(5)+B(2,2)*u1(6);
n1=[0;-1]; n4=[-1;0];

c(1)=n1'*B*(K1*v1+b1)+n1'*A*v1;
c(2)=n4'*B*(K1*v1+b1)+n4'*A*v1;
c(3)=n1'*xdot11;
c(4)=n4'*xdot21 ;
c(5)=K1*v1+b1-5;
c(6)=-K1*v1-b1+5;
ceq=[];

end

[tty]

yes K1 the 2x2 matrix and the vector b1 (1x2) should satisfy these constraints.

I tried to solve the system with x0=zeros(3,2) just to see if my prog will run and it does .But I want to solve the optimization at v1 wich has the coordinates in my state space and I don t know how to choose the x0?

clc;
close all;
v1=[80;60];
n1=[0;-1]; n4=[-1;0]; % normal vectors
A=[sqrt(2) 1;2 sqrt(3)];
B=[1 4;3 1];
nonlcon=@carreconst;
fun1=@(u1)(-(((A(1,1)+B(1,1)*u1(1)+B(1,2)*u1(3))*v1(1))+((A(1,2)+B(1,1)*u1(2)+B(1,2)*u1(4))*v1(2))+B(1, ...
  1)*u1(5)+B(1,2)*u1(6)));

x0=zeros(3,2);
sol=fmincon(fun1,x0,[],[],[],[],[],[],nonlcon);

m=sol((1:2),(1:2));
n=sol(3,:);
K1=m
 b1=n'
 tspa = [1 7];
iniCo = [80 60];
[t, x] = ode45(@(t,x) systLeft(t,x,K1,b1,A,B), tspa, iniCo');
plot(x(:,1),x(:,2)) 

function xdot = systLeft(t,x,KLeft,bLeft,A,B)

xdot = A*x + B*KLeft*x+B*bLeft;

end

[Walter Roberson]

Typically it is best to choose x0 as "something reasonable" for the solution. If you know that values near (say) 45 are common, then use those.

In general, it is often better not to have all of the elements of x0 be the same, and often better not to start at zero (or at least not all zero). In the case where the expressions happen to be of symmetric functions, using all elements the same can lead to problems of not giving a hint as to which direction to search; likewise it is common that using 0 as a start can be less.. informative... than a non-zero value.

[Torsten]

I must admit that I still don't understand

I want to find K1 and b1 such that I maximize xdot11 at a point v1 under these constraints :

What is v1 ? How does it correspond to something you get from or input into your differential equations ?

[tty]

my system should evolves inside a rectangle region in my state space and leave the rectangle only from the right edge (line between v2 and v3). this rectangle is specified with four vertices v1,v2,v3 and v4.

Answer 1

[Alan Weiss]

For an example that optimizes the solution of an ODE using optimization variables, see Fit ODE Parameters using Optimization Variables. That example shows calling ode45 to solve the ODE, and has an optimization problem built around this ODE solver.

For an entirely different approach, see Discretized Optimal Trajectory, Problem-Based. I don't know if this approach is relevant to your problem.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation