在 Matlab 中更新 fmincon 中的非线性约束数

Question

我正在尝试在while 循环中使用fmincon，这样在不满足while 条件之前，必须执行fmincon。每次fmincon 不能满足特定条件（例如x(N)-7.6==Tol）时，应更新（增加）非线性约束的数量N。 fmincon 怎么可能？ 假设我最初有 18 个非线性等式方程 (ceq(1) ... ceq(18))。当while条件不能满足时，下一次迭代时非线性等式方程的数量应增加到23个（ceq(1) ...ceq(23)）。

感谢您的创新想法....让我们为您提供有关我想做的事情的更多详细信息。我有一些非线性代数方程，所以我需要使用 NLP（非线性规划）求解器。除了我的成本函数是最小时间问题。实际上我的非线性约束方程是一些在时间坐标上离散化的动态控制方程。 N 是离散化的次数。基于拉格朗日优化技术，为了找到最优解，应将标量函数（拉格朗日）的梯度添加到系统中。正如我在我的问题中提到的，我需要用初始 N 来测试我的问题，然后如果 fmincon 的 Xoutput 不能满足约束，它需要增加描述的数量。它继续，直到 fmincon 的输出最佳答案足够接近我想要的答案。

Answer 1

剥开洋葱，在我看来，您的问题并没有解决您的真正问题，而且您正在以某种错误的方式从数字上解决您的问题。你如何以数字方式解决问题与你如何以分析方式解决问题（手工写出方程）相关但又不同。请在最后查看我的进一步评论，我建议与同学、实验室中的同学和/或教授交谈。

几个cmets：

1. 你的目标函数 @(x) norm(myfoctest(x)) 将总是返回 0，因为 myfoctest 返回空数组 [] 作为它的第一个参数，而在 Matlab 中，norm([]) 被定义为 0。

2. 您似乎打算解决minimize norm(f(x)) subject to f(x)==0 的问题，而不是minimize 0 subject to f(x)==0？我不明白在这种情况下约束 f(x)==0 的目的。为什么不minimize norm(f(x))？

3. 在你的函数中， 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的方式可能还有其他错误。

4. 其中许多约束是线性约束。将它们作为非线性约束输入会使 fmincon 的工作更加困难。

5. 我不知道最初的问题，但是这感觉你正在以一种随意、古怪的方式对它进行数值求解。我只看代码就发现了几个错误，如果我真的理解问题，我会担心还有更多错误。

6. 您有5*N+13 变量和5*N+13 非线性等式约束。你的可行集可能是一个点！题外话：很多优化算法都是从一个可行点开始，朝着可行的方向迈出一步。如果可行集是单点，则没有可行方向...在您的问题中，整个游戏都在寻找一个可行点（如果它甚至存在的话）？！

7. 我怀疑主要问题是您需要“在 Matlab 中更新 fmincon 中非线性约束的数量。”

8. 听起来好像您已经计算了拉格朗日的一阶条件，并且您正在输入这些作为优化问题的约束？如果是这样，那很可能不是你应该做的来解决数值问题。

一些建议...

1. 我认为您需要退后一步：以简洁的方式写下您要解决的问题，然后找出如何以有效的方式以数字方式解决它。我的直觉反应是，到目前为止，这有点混乱。

进一步评论

假设您有一些最小化问题（例如最优控制问题）：

minimize f(x) subject to g(x) <= 0.

其中 f 和 g 和凸，斯莱特条件成立，并且一阶条件是达到最小值的充分必要条件。你可以用数学方法解决这个问题并得到一些一阶条件：

dL/dx = 0

您可能认为用数值方法解决这个问题的方法是用数值方法求解方程组 dL/dx（来自 FOC）。如果 dL/dx 是一个线性方程组，这可能是正确的，但总的来说，这通常是解决问题的棘手方法。相反，您希望将f 和g 直接提供给您的优化算法。

要记住的一般要点：

• 求解线性方程组既高效又快速。
• 解决凸优化问题既高效又快速。
• 一般来说，求解非线性方程组或非凸优化可能是非常可怕的问题。

Answer 2

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

Answer 3

您有两种情况。首先，满足某些条件，其次，不满足。使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 将照常获得约束。但是如果条件不满足，约束条件就会改变。希望对您有所帮助。