解决最佳控制问题，ode45 vs fmincon

问题描述 投票：0回答：1

下午好！

首先，我寻找类似的问题一段时间，但（可能是因为我缺乏经验）我发现没有什么类似于我要问的问题。

我第一次使用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);

使用这段代码，我得到（显然）每个变量只有一个值作为结果，因为我没有找到任何涉及时间的方法，正如我之前所说，我开始寻找其他解决策略。

我发现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的基础教程;当然，这应该需要在各种文件中进行重构和排序（一个用于成本函数，一个用于约束）例如，但首先我想了解问题所在，然后我将以一个漂亮的形式。

非常感谢！

matlab compiler-errors ode45
1个回答
0
投票

通常你会把一些东西与Vector混淆在初始条件下，您声明了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的Symbolic Math Toolbox？

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.