解决最佳控制问题，ODE45与Fmincon

下午好！

首先，我一直在寻找类似的问题，但是（可能是因为我的经验不足）我没有发现与我要问的相似的问题。

我第一次使用MATLAB解决此类问题，因此我不确定该怎么做。简短的解释：

我正在为我的最佳控制课程做一个项目：我必须复制有关就业的论文的结果，并且我坚持了这些地块。我有以下数据：

五个变量函数（u（t），t（t），r（t），v1（t）和v2（t））
四个控制函数（U1（T），U2（T），U3（T），U4（T））
对控制变量的约束（每个U必须在0和1之间）
u，t，r，v1和v2的初始值（在t = 0中，尤其是v1和v2是恒定的）
哈密顿量λ系数的最终值

（注意：对于控件，我已经找到了最佳表达式，该表达式处于此形式：ui = min {1，max {0，" expression"}}。如果需要，我还可以给出四个表达，被忽略合成一点）

根据教授的建议，我试图使用fmincon，从理论上讲应该直接向我提供仅使用问题的成本函数来绘制一些结果的信息。但是在这种情况下，我有一些涉及计算时间的问题。以下是我用于FMINCON的代码：

syms u
%note: u(5) corresponds to U(t), but this is the only way I've found to get
%a result, the other u(i) are in ascending order (u(1) = u1 and so on...)
g = @(u) 30*u(5) + (20/2)*(u(1))^2 + (20/2)*(u(2))^2 + (10/2)*(u(3))^2 + (40/2)*(u(4))^2;
%initial guesses
u0 = [0 0 0 0 100000]; %
A = [];
b = [];
Aeq = [];
beq = [];
lb = 0.0 * ones(1,2,3,4);
ub = 1.0 * ones(1,2,3,4);
[x,fval,output,lambda] = fmincon(g, u0, A, b, Aeq, beq, lb, ub);

whit此代码，我（显然）每个变量仅获得一个值，并且由于我没有找到任何涉及时间的方法，所以我开始寻找其他求解策略。

我发现ODE45是一个微分方程求解器，它已包含在算法中的"时间迭代"，因此我尝试编写代码以使其与我的问题一起使用。我从论文中取出所有方程式，并将它们放入矢量中，如Mathworks示例所示，这是我的Matlab文件：

syms u1(t) u2(t) u3(t) u4(t)
syms U(t) T(t) R(t) V1(t) V2(t)
syms lambda_u lambda_t lambda_r lambda_v1 lambda_v2
%all the parameters provided by the paper
delta = 500;
alpha1 = 0.004;
alpha2 = 0.005;
alpha3 = 0.006;
gamma1 = 0.001;
gamma2 = 0.002;
phi1 = 0.22;
phi2 = 0.20;
delta1 = 0.09;
delta2 = 0.05;
k1 = 0.000003;
k2 = 0.000002;
k3 = 0.0000045;
%these two variable are set constant
V1 = 200;
V2 = 100;
%weight values for the cost function (only A1 is used in this case, but I left them all since the unused ones are irrelevant)
A1 = 30;
A2 = 20;
A3 = 20;
A4 = 10;
A5 = 40;
%ordering the unknowns in an array
x = [U T R u1 u2 u3 u4];
%initial conditions, ordered as the x vector (for the ui are guesses)
y0 = [100000 2000 1000 0 0 0 0];
%system set up
f = @(t,x) [delta - (1 + x(4))*k1*x(1)*V1 - (1 + x(5))*k2*x(1)*V2 - alpha1*x(1) + gamma1*x(2) + gamma2*x(3);...
    (1 + x(4))*k1*x(1)*V1 - k3*x(2)*V2 - alpha2*x(2) - gamma1*x(2);...
    (1 + x(5))*k2*x(1)*V2 - alpha3*x(3) - gamma2*x(3) + k3*x(2)*V2;...
    alpha2*x(2) + gamma1*x(2) + (1 + x(6))*phi1*x(1) + k3*x(2)*V2 - delta1*V1;...
    alpha3*x(3) + gamma2*x(3) + (1 + x(7))*phi2*x(1) - delta2*V2;...
    -A1 + (1 + x(4))*k1*V1*(lambda_u - lambda_t) + (1 + x(5))*k2*V2*(lambda_u - lambda_r) + lambda_u*alpha1 - lambda_v1*(1 + x(6))*phi1 - lambda_v2*(1 + x(7))*phi2;...
    -lambda_u*gamma1 + (alpha2 + gamma1)*(lambda_t - lambda_v1) + k3*V2*(lambda_t - lambda_r - lambda_v1);...
    -lambda_u*gamma2 + (alpha3 + gamma2)*(lambda_r - lambda_v2);...
    (1 + x(4))*k1*x(1)*(lambda_u - lambda_t) + lambda_v1*delta1;...
    (1 + x(5))*k2*x(1)*(lambda_u -lambda_r) + k3*x(2)*(lambda_t - lambda_r - lambda_v1) + lambda_v2*delta2];

%using ode45 to solve over the chosen time interval
[t,xa] = ode45(f,[0 10],y0);

使用此代码，我会收到以下错误：

Error using odearguments (line 95)
@(T,X)[DELTA-(1+X(4))*K1*X(1)*V1-(1+X(5))*K2*X(1)*V2-ALPHA1*X(1)+GAMMA1*X(2)+GAMMA2*X(3);(1+X(4))*K1*X(1)*V1-K3*X(2)*V2-ALPHA2*X(2)-GAMMA1*X(2);(1+X(5))*K2*X(1)*V2-ALPHA3*X(3)-GAMMA2*X(3)+K3*X(2)*V2;ALPHA2*X(2)+GAMMA1*X(2)+(1+X(6))*PHI1*X(1)+K3*X(2)*V2-DELTA1*V1;ALPHA3*X(3)+GAMMA2*X(3)+(1+X(7))*PHI2*X(1)-DELTA2*V2;-A1+(1+X(4))*K1*V1*(LAMBDA_U-LAMBDA_T)+(1+X(5))*K2*V2*(LAMBDA_U-LAMBDA_R)+LAMBDA_U*ALPHA1-LAMBDA_V1*(1+X(6))*PHI1-LAMBDA_V2*(1+X(7))*PHI2;-LAMBDA_U*GAMMA1+(ALPHA2+GAMMA1)*(LAMBDA_T-LAMBDA_V1)+K3*V2*(LAMBDA_T-LAMBDA_R-LAMBDA_V1);-LAMBDA_U*GAMMA2+(ALPHA3+GAMMA2)*(LAMBDA_R-LAMBDA_V2);(1+X(4))*K1*X(1)*(LAMBDA_U-LAMBDA_T)+LAMBDA_V1*DELTA1;(1+X(5))*K2*X(1)*(LAMBDA_U-LAMBDA_R)+K3*X(2)*(LAMBDA_T-LAMBDA_R-LAMBDA_V1)+LAMBDA_V2*DELTA2]
returns a vector of length 10, but the length of initial conditions vector is 7. The vector returned by
@(T,X)[DELTA-(1+X(4))*K1*X(1)*V1-(1+X(5))*K2*X(1)*V2-ALPHA1*X(1)+GAMMA1*X(2)+GAMMA2*X(3);(1+X(4))*K1*X(1)*V1-K3*X(2)*V2-ALPHA2*X(2)-GAMMA1*X(2);(1+X(5))*K2*X(1)*V2-ALPHA3*X(3)-GAMMA2*X(3)+K3*X(2)*V2;ALPHA2*X(2)+GAMMA1*X(2)+(1+X(6))*PHI1*X(1)+K3*X(2)*V2-DELTA1*V1;ALPHA3*X(3)+GAMMA2*X(3)+(1+X(7))*PHI2*X(1)-DELTA2*V2;-A1+(1+X(4))*K1*V1*(LAMBDA_U-LAMBDA_T)+(1+X(5))*K2*V2*(LAMBDA_U-LAMBDA_R)+LAMBDA_U*ALPHA1-LAMBDA_V1*(1+X(6))*PHI1-LAMBDA_V2*(1+X(7))*PHI2;-LAMBDA_U*GAMMA1+(ALPHA2+GAMMA1)*(LAMBDA_T-LAMBDA_V1)+K3*V2*(LAMBDA_T-LAMBDA_R-LAMBDA_V1);-LAMBDA_U*GAMMA2+(ALPHA3+GAMMA2)*(LAMBDA_R-LAMBDA_V2);(1+X(4))*K1*X(1)*(LAMBDA_U-LAMBDA_T)+LAMBDA_V1*DELTA1;(1+X(5))*K2*X(1)*(LAMBDA_U-LAMBDA_R)+K3*X(2)*(LAMBDA_T-LAMBDA_R-LAMBDA_V1)+LAMBDA_V2*DELTA2]
and the initial conditions vector must have the same number of elements.
Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in test (line 62)
[t,xa] = ode45(f,[0 10],y0);

我找不到解决方案，因为我使用了论文中给出的所有初始值。我剩下的唯一值是lambda系数的最终值，因为它们是最终值，我不确定是否可以使用它们。在这种情况下，我也不明白我应该将界限放在控制变量上。

有关完整性，我还将提供有关论文的链接：https://www.ripublication.com/ijss17/ijssv12n3_13.pdf

您能帮我弄清楚我能做什么来解决我的问题？

P.S：我知道这是一个很糟糕的代码，但我是基于Mathworks的基础教程；可以肯定的是，应该在各种文件中重新分配和订购这一点（一个用于成本函数，例如，一个用于约束），但首先我想了解问题所在，然后我将全部放置为漂亮的形式。<<<<<<<<<<<</p>

非常感谢！

通常，您将某些东西与向量相混淆。在初始条件下，您声明了7个值：

%initial conditions, ordered as the x vector (for the ui are guesses)
y0 = [100000 2000 1000 0 0 0 0];

，但您宣布了10个ODE：

%system set up
f = @(t,x) [delta - (1 + x(4))*k1*x(1)*V1 - (1 + x(5))*k2*x(1)*V2 - alpha1*x(1) + gamma1*x(2) + gamma2*x(3);...
    (1 + x(4))*k1*x(1)*V1 - k3*x(2)*V2 - alpha2*x(2) - gamma1*x(2);...
    (1 + x(5))*k2*x(1)*V2 - alpha3*x(3) - gamma2*x(3) + k3*x(2)*V2;...
    alpha2*x(2) + gamma1*x(2) + (1 + x(6))*phi1*x(1) + k3*x(2)*V2 - delta1*V1;...
    alpha3*x(3) + gamma2*x(3) + (1 + x(7))*phi2*x(1) - delta2*V2;...
    -A1 + (1 + x(4))*k1*V1*(lambda_u - lambda_t) + (1 + x(5))*k2*V2*(lambda_u - lambda_r) + lambda_u*alpha1 - lambda_v1*(1 + x(6))*phi1 - lambda_v2*(1 + x(7))*phi2;...
    -lambda_u*gamma1 + (alpha2 + gamma1)*(lambda_t - lambda_v1) + k3*V2*(lambda_t - lambda_r - lambda_v1);...
    -lambda_u*gamma2 + (alpha3 + gamma2)*(lambda_r - lambda_v2);...
    (1 + x(4))*k1*x(1)*(lambda_u - lambda_t) + lambda_v1*delta1;...
    (1 + x(5))*k2*x(1)*(lambda_u -lambda_r) + k3*x(2)*(lambda_t - lambda_r - lambda_v1) + lambda_v2*delta2];

上面代码中的每一行都被认为是一个ode。

但这不是全部。第二个问题是您的施工。您将符号数学（lambda称为syms）与数值求解混合在一起，这将很棘手。我不熟悉您要解决的确切科学问题，但是如果您无法避免象征性数学，也许您应该尝试dsolve中的CC_2？

相关内容

最新更新

热门标签：