我试图在while
循环中使用fmincon
,以便在不满足while
条件之前,必须执行fmincon
。每次fmincon
不能满足特定条件(例如,x(N)-7.6==Tol
(时,非线性约束的数量N
应该更新(增加(。这怎么可能fmincon
?假设我最初有 18 个非线性相等方程 ( ceq(1)
... ceq(18)
(。当不能满足while
条件时,非线性等式方程的数量应增加到 23 ( ceq(1)
... ceq(23)
( 在下一次迭代中。
TNX为您的创新理念...让我们为您提供有关我想做什么的更多详细信息。我有一组非线性代数方程,所以我需要使用NLP(非线性规划(求解器。除了我的成本函数是最小时间问题。实际上,我的非线性约束方程是一些在时间坐标中离散化的动态控制方程。N 是离散化的次数。基于拉格朗日优化技术,为了找到最优解,应在系统中增加标量函数梯度(拉格朗日函数(。正如我在问题中提到的,我需要用初始 N 测试我的问题,那么如果 fmincon 的 Xoutput 不能满足约束,它需要增加描述的数量。它继续,直到 fmincon 的输出最佳答案足够接近我想要的答案。
洋葱,在我看来,你的问题并没有解决你的真正问题,你正在以某种错误的方式用数字方式解决问题。你如何用数字解决问题是相关的,但与你如何用分析方式解决问题是相关的,但又是不同的,手写出方程式。请参阅我最后的进一步评论,我建议与同学、实验室的同学和/或教授交谈。
几点意见:
您的目标函数
@(x) norm(myfoctest(x))
将始终返回 0,因为 myfoctest 返回空数组[]
作为其第一个参数,而在 Matlab 中,norm([])
定义为 0。而不是
minimize 0 subject to f(x)==0
,你似乎打算解决问题minimize norm(f(x)) subject to f(x)==0
?在这种情况下,我不明白约束 f(x(==0 的目的。为什么不minimize norm(f(x))
?在您的函数中,
myfoctest
为什么Ceq(2*N)=x(5*N+12)-x(5*N+11)-x(4*N+2)
在 for 循环中?您将x(32)- x(31) - x(18)
的值分配给四个不同的时间Ceq(8)
(即 j=1:4(。这是你的意图吗?此错误向我表明,您编写myfoctest
的方式可能存在其他错误。其中许多约束是线性约束。将它们作为非线性约束输入将使 fmincon 的工作更加困难。
我不知道最初的问题,但在我看来,你正在以一种随意、糟糕的方式用数字方式解决它。我只是瞥了一眼代码就发现了几个错误,如果我真的理解了这个问题,我会担心会有更多错误。
您有
5*N+13
变量和5*N+13
非线性相等约束。你的可行集合可能是单点!题外话:许多优化算法从一个可行的点开始,朝着可行的方向迈出一步。如果可行集合是单个点,则没有可行的方向...在你的问题中,整个游戏都在找到一个可行的点(如果它甚至存在(?!我怀疑主要问题是你需要"更新 Matlab 中 fmincon 中的非线性约束数量"。
听起来好像您已经计算了拉格朗日量的一阶条件,并且您正在输入这些条件作为优化问题的约束?如果是这样,那很可能不是您应该做的数值解决。
一些建议...
- 我认为你需要回到起点:以干净的方式写下你试图解决的问题,然后弄清楚如何以有效的方式以数字方式解决它。我的直觉反应是,到目前为止,这是一个有点混乱的混乱。
进一步评论
假设您有一些最小化问题(例如,最佳控制问题(:
minimize f(x) subject to g(x) <= 0.
其中 f 和 g 和凸,斯莱特条件成立,一阶条件是达到最小值的必要条件和充分条件。你可以用数学方法解决这个问题,得到一些一阶条件:
dL/dx = 0
您可能会认为以数值方式解决此问题的方法是以数值方式求解方程组 dL/dx(来自 FOC(。如果dL/dx是一个线性方程组,这可能是真的,但总的来说,这通常是解决问题的棘手方法。相反,您希望将f
和g
直接馈送到优化算法中。
要记住的一般要点:
- 求解线性方程组既高效又快速。
- 求解凸优化问题既高效又快速。
- 一般来说,求解非线性方程组或非凸优化可能是可怕的问题。
你有两种情况。首先,满足某些条件,其次,则不满足。使while
语句在每种情况下都成立。在循环中添加标志变量,如果不满足条件,该变量将更改其值。例如,您可以放置以下内容:
flag = (x(N)-7.6<Tol);
如果满足条件,这将返回1
,否则0
。
在mycon
函数中,将flag
添加为输入变量:
function [c,ceq] = mycon(all_variables_you_had_before,flag)
然后,添加mycon
中的逻辑块,如下所示:
if flag == 1
ceq = [___]; %//put your 18 conditions here
else
ceq = [___]; %//put your 23 conditions here
end
最后,不要忘记在主脚本的fmincon
行中添加mycon(all_variables_you_had_before,flag)
:
x = fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@(all_variables_you_had_before) mycon(all_variables_you_had_before,flag))
因此,如果满足条件,您的fmincon
将像往常一样获得约束。但是,如果不满足条件,约束就会改变。希望有帮助。
function [x]=runnested(x0,N)
r=ones(4,1);
N=length(r);
Tol=0.001;
for k=1:N
for i=1:N
x0=rand(5*N+13,1)
options = optimset('Largescale','off','algorithm','interior-point','Display','iter');
[x(i,:),fval,exitflag,output]=fmincon(@(x) norm(myfoctest(x)),x0,[],[],[],[],[],[],@myfoctest,options)
end
if x(N)-7.61<=Tol
break;
else
N=N+1;
end
end
function [C,Ceq]=myfoctest(x,N,r)
C=[];
r=ones(4,1);
N=length(r);
f=3.5e-6; %km/s^2
i1=10*(pi/180);
Ts=110; %sec
V0=7.79; %km/sec
a1=7.61; %km/sec
b1=0.01*a1;
a2=20*(pi/180); % rad %10 deg
b2=0.01*a2; %rad
Omeg0=10*(pi/180); %rad
Ceq=zeros(5*N+13,1);
for j=1:N-1
Ceq(j)=x(3*N+1+j)- x(3*N+j)-2*x(4*N+1+j)*Ts*f*sin(x(2*N+1+j))./(pi*sin(i1)*x(j)^2)
Ceq(N)=x(5*N+10)-x(5*N+9)-x(3*N+2) %x(5*N+10)-x(5*N+9)-x(4*N+7)
Ceq(N+j)=x(4*N+1+j)-x(4*N+j)
Ceq(2*N)=x(5*N+12)-x(5*N+11)-x(4*N+2)
Ceq(2*N+1)=x(3*N+1)*Ts*f*sin(x(2*N+1))+2*x(4*N+1)*Ts*f*cos(x(2*N+1))/(pi*V0*sin(i1))
Ceq(2*N+1+j)=x(3*N+1+j)*Ts*f*sin(x(2*N+1+j))+2*x(4*N+1+j)*Ts*f*cos(x(2*N+1+j))./(pi*x(j)*sin(i1))
Ceq(3*N+1)=1-x(5*N+9)*b1-x(5*N+10)*b1-x(5*N+11)*b2-x(5*N+12)*b2-x(5*N+8)*N*Ts/100-x(5*N+13)
Ceq(3*N+2)=-2*x(5*N+8)*x(5*N+2)
Ceq(3*N+3)=-2*x(5*N+9)*x(5*N+3)
Ceq(3*N+4)=-2*x(5*N+10)*x(5*N+4)
Ceq(3*N+5)=-2*x(5*N+11)*x(5*N+5)
Ceq(3*N+6)=-2*x(5*N+12)*x(5*N+6)
Ceq(3*N+7)=2*x(5*N+13)*cos(x(5*N+7))*sin(x(5*N+7))
Ceq(3*N+8)=V0-x(1)-Ts*f*cos(x(2*N+1))
Ceq(3*N+8+j)=x(j)-x(j+1)-Ts*f*cos(x(2*N+1+j))
Ceq(4*N+8)=Omeg0-x(N+1)+2*Ts*f*sin(x(2*N+1))/(pi*V0*sin(i1))
Ceq(4*N+8+j)=Omeg0-x(j+1)+2*Ts*f*sin(x(2*N+1+j))./(pi*x(j)*sin(i1))
Ceq(5*N+8)=-x(5*N+2)^2-N*Ts/100-N*Ts*x(3*N+1)/100
Ceq(5*N+9)=-x(5*N+3)^2-x(N)+a1+b1-b1*x(3*N+1)+7.61/100
Ceq(5*N+10)=-x(5*N+4)^2+x(N)+a1+b1-b1*x(3*N+1)-7.61/100
Ceq(5*N+11)=-x(5*N+5)^2-x(2*N)+a2+b2-b2*x(3*N+1)+0.35/100
Ceq(5*N+12)=-x(5*N+6)^2+x(2*N)+a2+b2-b2*x(3*N+1)-0.35/100
Ceq(5*N+13)=-(sin(x(5*N+7)))^2-x(5*N+1)
end
end
end