Minimize a variable in a no-linear system

2023 1 21
1 个回答
回答已采纳
更新时间：2023 1 31
5 次查看（30 天）

0 个投票

Hi, i need to minimize Fv in this system:
The problems variables are :

19 个评论
显示 17更早的评论 隐藏 17更早的评论

David Hill
David Hill 2023-1-21
Show us your code and ask a specific question.
Torsten
Torsten 2023-1-21
编辑：Torsten 2023-1-21
Use "fmincon".
Matt J
Matt J 2023-1-22
I would imagine linprog would be applicable, since all constraints appear to be linear.
Torsten
Torsten 2023-1-22
Except the second to last.
Antonio Cassano
Antonio Cassano 2023-1-22
编辑：Matt J 2023-1-23
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
Torsten 2023-1-22
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
Antonio Cassano 2023-1-23
编辑：Antonio Cassano 2023-1-23
So isn't there a method to solve this problem in Matlab?
Antonio Cassano
Antonio Cassano 2023-1-23
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
Matt J 2023-1-23
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
Matt J 2023-1-23
编辑：Matt J 2023-1-23
Incidentally, you could have had,
problem.Objective = Fv;
Antonio Cassano
Antonio Cassano 2023-1-26
编辑：Matt J 2023-1-26
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
Matt J 2023-1-26
Again, you appear to have more equaltity constraints than unknowns, so the problem is extremely likely to be infeasible.
Antonio Cassano
Antonio Cassano 2023-1-29
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
Torsten 2023-1-29
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
Antonio Cassano 2023-1-30
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
Matt J 2023-1-30
Well, maybe it eventually will converge.
Antonio Cassano
Antonio Cassano 2023-1-31
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
Matt J 2023-1-31
@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
Matt J 2023-1-31
编辑：Matt J 2023-1-31
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.

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Matt J
Matt J 2023-1-22
编辑：Matt J 2023-1-22

0 个投票

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...

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

产品

版本

R2022b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by