Minimize a variable in a no-linear system

Hi, i need to minimize Fv in this system:

The problems variables are :

[David Hill]

Show us your code and ask a specific question.

[Torsten]

Use "fmincon".

[Matt J]

I would imagine linprog would be applicable, since all constraints appear to be linear.

[Torsten]

Except the second to last.

[Antonio Cassano]

Re=0.055;
densita=75;
volume_oggetto=0.091875;
G=[0.117 0 -0.292];
beta=0:90;
n=length(beta);
g=9.807;
pa=101325;
fst=[0.4 , 0.3 , 0.25];
Fvt=zeros(3,n);
p0t=zeros(3,n);
x_kt=zeros(3,n);
y_kt=zeros(3,n);

for i=1:n
    for j=1:1
        angolo=deg2rad(i-1);
        
        % Create optimization variables
        Fv = optimvar("Fv",'LowerBound',0,'UpperBound',962);
        Fn1 = optimvar("Fn1",'LowerBound',0,'UpperBound',10000);
        Fn2 = optimvar("Fn2",'LowerBound',0,'UpperBound',10000);
        Fn3 = optimvar("Fn3",'LowerBound',0,'UpperBound',10000);
        Fn4 = optimvar("Fn4",'LowerBound',0,'UpperBound',10000);
        Ftx1 = optimvar("Ftx1",'LowerBound',-2000,'UpperBound',2000);
        Ftx2 = optimvar("Ftx2",'LowerBound',-2000,'UpperBound',2000);
        Ftx3 = optimvar("Ftx3",'LowerBound',-2000,'UpperBound',2000);
        Ftx4 = optimvar("Ftx4",'LowerBound',-2000,'UpperBound',2000);
        Fty1 = optimvar("Fty1",'LowerBound',-2000,'UpperBound',2000);
        Fty2 = optimvar("Fty2",'LowerBound',-2000,'UpperBound',2000);
        Fty3 = optimvar("Fty3",'LowerBound',-2000,'UpperBound',2000);
        Fty4 = optimvar("Fty4",'LowerBound',-2000,'UpperBound',2000);
        x_k = optimvar("x_k",'LowerBound',-2,'UpperBound',2);
        y_k = optimvar("y_k",'LowerBound',-2,'UpperBound',2);
        % Set initial starting point for the solver
        initialPoint.Fv = repmat(300,size(Fv));
        initialPoint.Fn1 =  zeros(size(Fn1));
        initialPoint.Fn2 =  zeros(size(Fn2));
        initialPoint.Fn3 =  zeros(size(Fn3));
        initialPoint.Fn4 = zeros(size(Fn4));
        initialPoint.Ftx1= repmat(1,size(Ftx1));
        initialPoint.Ftx2 = repmat(1,size(Ftx2));
        initialPoint.Ftx3 = repmat(1,size(Ftx3));
        initialPoint.Ftx4 = repmat(1,size(Ftx4));
        initialPoint.Fty1 = zeros(size(Fty1));
        initialPoint.Fty2 = zeros(size(Fty2));
        initialPoint.Fty3 = zeros(size(Fty3));
        initialPoint.Fty4 = zeros(size(Fty4));
        initialPoint.x_k = zeros(size(x_k));
        initialPoint.y_k = repmat(0.01,size(y_k));
        % Create problem
        problem = optimproblem;
        % Define problem objective
        problem.Objective = fcn2optimexpr(@forzadivuoto,Fv);
        % Define problem constraints
        problem.Constraints.constraint1 = Fv-Fn1 -Fn2-Fn3-Fn4-densita*g*volume_oggetto*cos(angolo) == 0;
        problem.Constraints.constraint2 = Ftx1+Ftx2+Ftx3+Ftx4==0 ;
        problem.Constraints.constraint3 = -densita*volume_oggetto*g*sin(angolo)+Fty1+Fty2+Fty3+Fty4 == 0;
        problem.Constraints.constraint4 = -densita*volume_oggetto*g*sin(angolo)*abs(G(1,3))-Re*(Fn4-Fn2) == 0;
        problem.Constraints.constraint5 = -densita*volume_oggetto*g*cos(angolo)*abs(G(1,1))-Re*(Fn3-Fn1) == 0;
        problem.Constraints.constraint6 = -densita*volume_oggetto*g*G(1,1)*sin(angolo)*-Re*(Ftx4-Ftx2+Fty3-Fty1) == 0 ;
        problem.Constraints.constraint7 = Fn1+Fn3==Fn2+Fn4 ;
        problem.Constraints.constraint8 = Fty1/Ftx1==-x_k/y_k ;
        problem.Constraints.constraint9 = Fty2/Ftx2==-x_k/y_k ;
        problem.Constraints.constraint10 = Fty3/Ftx3==-x_k/y_k ;
        problem.Constraints.constraint11 = Fty4/Ftx4==-x_k/y_k ;
        problem.Constraints.constraint12 = Ftx1^2+Fty1^2==fst(1,j)*Fn1^2 ;
        problem.Constraints.constraint13 = Ftx2^2+Fty2^2==fst(1,j)*Fn2^2 ;
        problem.Constraints.constraint14 = Ftx3^2+Fty3^2==fst(1,j)*Fn3^2 ;
        problem.Constraints.constraint15 = Ftx4^2+Fty4^2==fst(1,j)*Fn4^2 ;
        % problem.Constraints.constraint16 = Fn1 >= 0;
        % problem.Constraints.constraint17 = Fn2 >= 0;
        % problem.Constraints.constraint18 = Fn3 >= 0;
        % problem.Constraints.constraint19 = Fn4 >= 0;
        
        % Set nondefault solver options
        
        options = optimoptions("fmincon","PlotFcn",["optimplotfval","optimplotconstrviolation","optimplotfirstorderopt","optimplotstepsize","optimplotfvalconstr","optimplotx"],'Algorithm','interior-point');
        
        % Display problem information
        
        % show(problem);
        
        % Solve problem
        [solution,objectiveValue,reasonSolverStopped] = solve(problem,initialPoint,"Solver","fmincon","Options",options);
        % Display results
        disp(solution)
        disp(reasonSolverStopped)
        disp(objectiveValue)
        % Clear variables
        clearvars Fv Fn1 Fn2 Fn3 Fn4 Ftx1 Ftx2 Ftx3 Ftx4 Fty1 Fty2 Fty3 Fty4 x_k y_k initialPoint options
        
        Fvt(j,i) = solution.Fv ;
        p0t(j,i) = solution.Fn1 ;
        x_kt(j,i) = solution.x_k ;
        y_kt(j,i) = solution.y_k ;
    end
end

% Objective function 
function objective = forzadivuoto(Fv)
objective = Fv;
end

This is my code that give me only infeasible point.

[Torsten]

If I counted correctly, you have 15 unknowns and 15 equality constraints.

So you can be happy if you even get a solution for your problem. Bound constraints on the solution can not be imposed.

[Antonio Cassano]

So isn't there a method to solve this problem in Matlab?

[Antonio Cassano]

This is a semplification of my original problem where i have only 4 variables but there are some difficult integrals to examine,if you want i can send my original script for the real problem

[Matt J]

There is nothing more to say about how to solve the problem. The only question is, does the problem have a solution? The 15-variable fake problem apparently does not have a solution, but it might not matter in the end, because maybe your real problem does.

[Matt J]

Incidentally, you could have had,

problem.Objective = Fv;

[Antonio Cassano]

This is the script for my real problem but this don't give me a faesible solution.

Re=0.055;
densita=75;
volume_oggetto=0.091875;
G=[0.117 0 -0.292];
beta=0:90;
n=length(beta);
g=9.807;
pa=101325;
fst=[0.4 , 0.3 , 0.25];


Tx=zeros(1,n);
Ty=zeros(1,n);
Tz=zeros(1,n);
a=zeros(1,n);
b=zeros(1,n);

Fvt=zeros(3,n);
p0t=zeros(3,n);
x_kt=zeros(3,n);
y_kt=zeros(3,n);
ft=zeros(3,n);




%%
for i=1:n

    angolo=deg2rad(i-1);
    
    
    Tx(i)=-densita*volume_oggetto*g*sin(angolo)*abs(G(1,3));
    Ty(i)=densita*volume_oggetto*g*cos(angolo)*G(1,1);
    Tz(i)=-densita*volume_oggetto*g*sin(angolo)*G(1,1);
    a(i)=-Ty(i)/(pi*Re^3);
    b(i)=Tx(i)/(pi*Re^3);
    
end



for i=1:n

    for j=1:1

    syms f1 x_k1 y_k1 p01 x

    angolo=deg2rad(i-1);

    fun1=f1*Re*(p01+Re*(a(i)*cos(x)+b(i)*sin(x)))*(y_k1-Re*sin(x))/sqrt((x_k1-Re*cos(x))^2+(y_k1-Re*sin(x))^2)  ;

     f1fnc=matlabFunction(fun1,'var',{x,f1,p01,x_k1,y_k1});
     F1=int(f1fnc,x,0,2*pi);
     F1fnc_=matlabFunction(F1,'var',{f1,p01,x_k1,y_k1});

    fun2=-f1*Re*(p01+Re*(a(i)*cos(x)+b(i)*sin(x)))*(x_k1-Re*cos(x))/sqrt((x_k1-Re*cos(x))^2+(y_k1-Re*sin(x))^2) ;
 f2fnc=matlabFunction(fun2,'var',{x,f1,p01,x_k1,y_k1});
     F2=int(f2fnc,x,0,2*pi);
     F2fnc_=matlabFunction(F2,'var',{f1,p01,x_k1,y_k1});

    fun3=Re^2*(-f1*(p01+Re*(a(i)*cos(x)+b(i)*sin(x)))*cos(x)*(x_k1-Re*cos(x))/sqrt((x_k1-Re*cos(x))^2+(y_k1-Re*sin(x))^2)-f1*(p01+Re*(a(i)*cos(x)+b(i)*sin(x)))*sin(x)*(y_k1-Re*sin(x))/sqrt((x_k1-Re*cos(x))^2+(y_k1-Re*sin(x))^2));
      f3fnc=matlabFunction(fun3,'var',{x,f1,p01,x_k1,y_k1});
     F3=int(f3fnc,x,0,2*pi);
     F3fnc_=matlabFunction(F3,'var',{f1,p01,x_k1,y_k1});

     

% Create optimization variables
Fv1 = optimvar("Fv");
p01 = optimvar("p0");
f1 = optimvar("f");
x_k1 = optimvar("x_k");
y_k1 = optimvar("y_k");

% Set initial starting point for the solver

if i==1

initialPoint.Fv = repmat(350,size(Fv1));
initialPoint.p0 = repmat(1500,size(p01));
initialPoint.f = repmat(0.0001,size(f1));
initialPoint.x_k = repmat(-0.0001,size(x_k1));
initialPoint.y_k = repmat(-0.0001,size(y_k1));

else 
    
initialPoint.Fv = Fvt(j,i-1);
initialPoint.p0 = p0t(j,i-1);
initialPoint.x_k = x_kt(j,i-1);
initialPoint.y_k = y_kt(j,i-1);
initialPoint.f = ft(j,i-1);

end




% Create problem
problem = optimproblem;

% Define problem objective
problem.Objective = Fv1;

F1eq=fcn2optimexpr(F1fnc_,f1,p01,x_k1,y_k1);
F2eq=fcn2optimexpr(F2fnc_,f1,p01,x_k1,y_k1);
F3eq=fcn2optimexpr(F3fnc_,f1,p01,x_k1,y_k1);

% Define problem constraints

problem.Constraints.constraint1 = p01 >= 0;
problem.Constraints.constraint2 = F1eq == 0;
problem.Constraints.constraint3 = -densita*volume_oggetto*g*sin(angolo)+F2eq == 0;
problem.Constraints.constraint4 = Fv1 - 2*p01*pi*Re - densita*g*volume_oggetto*cos(angolo) == 0;
problem.Constraints.constraint5 = f1 <= fst(1,j);
problem.Constraints.constraint6 = F1eq.^2+F2eq.^2-f1.^2*(2*pi*p01*Re).^2 == 0;
problem.Constraints.constraint7 = F2eq ./ F1eq == -x_k1/y_k1;
problem.Constraints.constraint8 = -densita*volume_oggetto*g*G(1,1)*sin(angolo)+F3eq == 0 ;

% Set nondefault solver options

options = optimoptions("fmincon","PlotFcn",["optimplotfval","optimplotconstrviolation","optimplotfirstorderopt","optimplotstepsize","optimplotfvalconstr","optimplotx"]);

% Display problem information

show(problem);

% Solve problem

[solution,objectiveValue,reasonSolverStopped] = solve(problem,initialPoint,"Solver","fmincon","Options",options);

% Display results

disp(solution)
disp(reasonSolverStopped)
disp(objectiveValue)


 Fvt(j,i) = solution.Fv ;
 p0t(j,i) = solution.p0 ;
 x_kt(j,i) = solution.x_k ;
 y_kt(j,i) = solution.y_k ;
 ft(j,i)=solution.f;


% Clear variables

clearvars Fv1 p01 x_k1 y_k1 f1 initialPoint options
     
    end
end

[Matt J]

Again, you appear to have more equaltity constraints than unknowns, so the problem is extremely likely to be infeasible.

[Antonio Cassano]

Sorry but if i have n unknowns with n-1 equations, with an optimization algorithm i can find the optimal solution because the algorithm hypothesizes the optimization variable therefore it's possible to solve the system of n-1 unk and n-1 eq.

Right?

[Torsten]

5 unknows:

initialPoint.Fv = Fvt(j,i-1);
initialPoint.p0 = p0t(j,i-1);
initialPoint.x_k = x_kt(j,i-1);
initialPoint.y_k = y_kt(j,i-1);
initialPoint.f = ft(j,i-1);

6 equality constraints, 2 inequality constraints:

problem.Constraints.constraint1 = p01 >= 0;
problem.Constraints.constraint2 = F1eq == 0;
problem.Constraints.constraint3 = -densita*volume_oggetto*g*sin(angolo)+F2eq == 0;
problem.Constraints.constraint4 = Fv1 - 2*p01*pi*Re - densita*g*volume_oggetto*cos(angolo) == 0;
problem.Constraints.constraint5 = f1 <= fst(1,j);
problem.Constraints.constraint6 = F1eq.^2+F2eq.^2-f1.^2*(2*pi*p01*Re).^2 == 0;
problem.Constraints.constraint7 = F2eq ./ F1eq == -x_k1/y_k1;
problem.Constraints.constraint8 = -densita*volume_oggetto*g*G(1,1)*sin(angolo)+F3eq == 0 ;

Or did I interprete something wrong ?

[Antonio Cassano]

Yes, but i can remove constraint6 and 7 therefore i will have 5 unks and 4 equation (remain only equilibrium equation);

i tried this but the algorithm didn't converge yet.

[Matt J]

Well, maybe it eventually will converge.

[Antonio Cassano]

Well, the algorithm converge but in the initial point the FirstOrderOptm is 0.3 therefore i have a jump of my variable, how can i improve it?

[Matt J]

@Antonio Cassano I would post that as a new question. First of all, we don't know what your current code looks like, the one that produced that graph. Secondly, we don't know what you would consider an "improvement" to the current result. All that could and should go into a new post..

[Matt J]

However you might also try looping backwards from n:-1:1 and initializing instead with,

initialPoint.Fv = Fvt(j,i+1);
initialPoint.p0 = p0t(j,i+1);
initialPoint.x_k = x_kt(j,i+1);
initialPoint.y_k = y_kt(j,i+1);
initialPoint.f = ft(j,i+1);

since the solutions seem to be more sensitive to the initial guess at lower i.

Answer 1

[Matt J]

If you use the Problem Based optimization tools, you can pretty much type in the problem as you have for us above, and the solver will be selected for you.

Fv=optimvar('Fv',4,1);
Fn=optimvar('Fn',4,1);
Ftx=optimvar('Ftx',4,1);
Fty=optimvar('Fty',4,1);


con(1) = Fez+4*Fv-sum(Fn)==0;

con(2) = Fex+sum(Ftx)==0;

etc...

[Antonio Cassano]

I tried this solution but the solver give me only infeasible point.