matlab - 解决最佳控制问题，ode45与fmincon

Question

下午好!

首先，我找了一段时间类似的问题，但是(可能是由于我的经验不足)我没有发现与我要问的问题类似的东西。

我是第一次使用matlab解决此类问题，因此我不确定该怎么做。简要说明:

我正在为我的“最佳控制”类(class)做一个项目:我必须复制一份关于就业的论文的结果，而我一直陷于困境。我有以下数据:

五个变量函数(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是一个微分方程求解器，其算法中已经包含了“时间迭代”，因此我尝试编写代码以使其能够解决我的问题。
我从论文中提取了所有方程式，并将其放入 vector 中，如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

您能帮我弄清楚该如何解决我的问题吗？

附注:我知道这是一个非常糟糕的代码，但我基于mathworks的基础教程；确保应该将其重构并排序到各种文件中(例如，一个用于成本函数，一个用于约束)，但是首先，我想了解问题出在哪里，然后将其以漂亮的形式显示。

非常感谢!

Answer 1

通常，您会将 vector 与 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吗？