You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Minimize a variable in a no-linear system
1 view (last 30 days)
Show older comments
Antonio Cassano
on 21 Jan 2023
Hi, i need to minimize Fv in this system:
The problems variables are :
19 Comments
David Hill
on 21 Jan 2023
Show us your code and ask a specific question.
Antonio Cassano
on 22 Jan 2023
Edited: Matt J
on 23 Jan 2023
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.
Antonio Cassano
on 23 Jan 2023
Edited: Antonio Cassano
on 23 Jan 2023
So isn't there a method to solve this problem in Matlab?
Antonio Cassano
on 23 Jan 2023
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
on 23 Jan 2023
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.
Antonio Cassano
on 26 Jan 2023
Edited: Matt J
on 26 Jan 2023
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
Antonio Cassano
on 29 Jan 2023
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
on 29 Jan 2023
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
on 30 Jan 2023
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.
Antonio Cassano
on 31 Jan 2023
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
on 31 Jan 2023
@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
on 31 Jan 2023
Edited: Matt J
on 31 Jan 2023
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.
Accepted Answer
Matt J
on 22 Jan 2023
Edited: Matt J
on 22 Jan 2023
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...
1 Comment
Antonio Cassano
on 22 Jan 2023
I tried this solution but the solver give me only infeasible point.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)