0

下午好！

首先，我找了一段时间类似的问题，但是（可能是因为我的经验不足）我没有发现任何与我要问的类似的问题。

我是第一次使用matlab来解决这类问题，所以不知道该怎么办。简要说明：

我正在为我的 Optimal Control 课程做一个项目：我必须复制一篇关于就业的论文的结果，但我被情节困住了。我有以下数据：

  • 五个变量函数（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

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

PS：我知道这是一个非常糟糕的代码，但我基于 mathworks 的基础教程；当然，这应该需要在各种文件中进行重构和排序（例如，一个用于成本函数，一个用于约束），但首先我想了解问题出在哪里，然后我会将所有文件都放在一个漂亮的形式中。

太感谢了！

4

1 回答 1

0

通常，您会将某些内容与向量混淆。在初始条件下，您声明了 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?

于 2019-03-15T18:54:02.183 回答