用MATLAB解《数理战术学》中的微分方程

这篇文章摘抄自沙昌基教授《数理战术学》第七章 多兵种对多兵种作战的微分对策模型，提供了求解微分方程的代码，防止自己忘记。

【例】最优火力分配矩阵随时间变化的例子

设一对二作战中双方的实力向量为

$$x=[x_1],y=\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]$$
其初值为
$$x_{10}=100,y_{10}=30,y_{20}=30$$
双方毁伤系数矩阵为
$$A=\left[\begin{array}{c}1\\ 1\end{array}\right],B=\left[\begin{array}{cc}9& 1\end{array}\right]$$
火力分配矩阵为
$$\Phi =\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right],\Psi =\left[\begin{array}{cc}{\psi }_1& {\psi }_2\end{array}\right]$$
其中
$$\begin{gathered} 0\le {\psi }_1\le 1,0\le {\psi }_2\le 1, \\ 0\le {\varphi }_1,0\le {\varphi }_2,{\varphi }_1+{\varphi }_2\le 1 \end{gathered}$$
于是，双方交战的Lanchester方程变为
$$\{\begin{array}{c}\begin{array}{rl}{\dot{x}}_1& =-9{\psi }_1{y}_1-{\psi }_2{y}_2,\\ {\dot{y}}_1& =-{\varphi }_1{x}_1,\\ {\dot{y}}_2& =-{\varphi }_2{x}_1.\end{array}\end{array}\text{\hspace{1em}(1)}$$
设双方的作战指数向量为
$$u=(9),v=\left[\begin{array}{c}1\\ 9\end{array}\right]$$
战斗时间为
$$T=\frac{1}{6}\mathrm{l}\mathrm{n}10+\mathrm{l}\mathrm{n}\frac{10}{9}$$
则支付函数为
$$\Theta =\Theta (x_1(T),y_1(T),y_2(T))=9x_1(T)-y_1(T)-9y_2(T)\text{\hspace{1em}(2)}$$

---

【解】

（1）构建Hamilton函数，为
$$H=\lambda (-9{\psi }_1{y}_1-{\psi }_2{y}_2)+{\mu }_1(-{\varphi }_1{x}_1)+{\mu }_2(-{\varphi }_2{x}_1)\text{\hspace{1em}(3)}$$
（2）根据Hamilton函数，求得协态方程为
$$\{\begin{array}{c}\begin{array}{rl}\dot{\lambda }& =-\frac{\partial H}{\partial x}={\mu }_1{\varphi }_1+{\mu }_2{\varphi }_2,\\ {\dot{\mu }}_1& =-\frac{\partial H}{\partial y_1}=9\lambda {\psi }_1,\\ {\dot{\mu }}_2& =-\frac{\partial H}{\partial y_2}=\lambda {\psi }_2\end{array}\end{array}\text{\hspace{1em}(4)}$$
（3）横截条件为
$$\{\begin{array}{c}\begin{array}{rl}\lambda (T)& =\frac{\partial \Theta }{\partial x_1}{|}_T=9,\\ {\mu }_1(T)& =\frac{\partial \Theta }{\partial y_1}{|}_T=-1,\\ {\mu }_2(T)& =\frac{\partial \Theta }{\partial y_2}{|}_T=-9\end{array}\end{array}\text{\hspace{1em}(5)}$$
参考单兵种对多兵种作战的微分对策模型对于协态变量的分析可知，恒有
$$\lambda (t)>0,{\mu }_1(t)<0,{\mu }_2(t)<0,0\le t\le T.$$
由鞍点条件可知，对于最优火力分配矩阵有
$$\{\begin{array}{c}{\psi }_1(t)={\psi }_2(t)=1,0\le t\le T,\\ {\varphi }_1(t)+{\varphi }_1(t)=1,0\le t\le T\end{array}\text{\hspace{1em}(6)}$$
且当$-{\mu }_1(t)>-{\mu }_2(t)$时，有
$${\varphi }_1(t)=1,{\varphi }_2(t)=0.\text{\hspace{1em}(7)}$$
当${\mu }_1(t)<-{\mu }_2(t)$时，有
$${\varphi }_1(t)=0,{\varphi }_2(t)=1.\text{\hspace{1em}(8)}$$
今因$-{\mu }_1(T)<-{\mu }_2(T)$，故在临近时刻 $T$ 的一个时间区间内上式均成立，即存在$\Delta \in [0,T)$，使得
$$-{\mu }_1(t)<-{\mu }_2(t),\Delta <t\le T\text{\hspace{1em}(9)}$$
现在来考察$\Delta$的最小值，并将这最小值仍记作$\Delta$。如果该$\Delta >0$，则意味着$\Delta$再不能减小，从而必有
$$-{\mu }_1(\Delta )=-{\mu }_2(\Delta )\text{\hspace{1em}(9)}$$
由式$(8)$有
$${\varphi }_1(t)=0,{\varphi }_2(t)=1,\Delta <t\le T$$
那么，联立式$(4)$、式$(6)$和式$(8)$，则式$(4)$变为
$$\{\begin{array}{c}\begin{array}{rl}\dot{\lambda }& ={\mu }_2,\\ {\dot{\mu }}_1& =9\lambda ,\\ {\dot{\mu }}_2& =\lambda ,\end{array}\end{array}\Delta <t\le T$$
考虑终值条件式$(5)$进行求解。

此处编写MATLAB代码求解上述微分方程，仿真代码为

%% 《数理战术学》第7.3节 式（7.23）、式（7.29）
% 【delta < t <= T】 区间
syms lambda(t) mu1(t) mu2(t) T
eq1 = diff(lambda,t)==mu2;con1 = lambda(T)==9;
eq2 = diff(mu1,t)==9*lambda;con2 = mu1(T)==-1;
eq3 = diff(mu2,t)==lambda;con3 = mu2(T)==-9;
sol1 = dsolve(eq1,eq2,eq3,con1,con2,con3);

求解结果为

书中求解结果为
$$\{\begin{array}{c}\begin{array}{rl}\lambda & =9e^{T-t},\\ {\mu }_1& =-81e^{T-t}+80,\\ {\mu }_2& =-9e^{T-t},\end{array}\end{array}\text{\hspace{1em}(10)}$$
MATLAB求解结果与书中求解结果相同。

下面来求$\Delta$，因式$(10)$在$\Delta \le t\le T$中成立，且$\Delta$不能再小，故由式$(9)$可得
$$\begin{array}{rl}-81e^{T-\Delta }+80& =-9e^{T-\Delta }\\ \Rightarrow e^{T-\Delta }& =\frac{10}{9}\\ \Rightarrow T-\Delta & =\mathrm{l}\mathrm{n}\frac{10}{9}\\ \Rightarrow \Delta & =T-\mathrm{l}\mathrm{n}\frac{10}{9}\\ & =\frac{1}{6}\mathrm{l}\mathrm{n}10\end{array}\text{\hspace{1em}(11)}$$
此处求解$\Delta$的MATLAB代码为

%% 《数理战术学》第7.3节 式（7.30）
syms delta
T = (1/6)*log(10) + log((10/9));
res.lambda = subs(sol1.lambda, t, delta);
res.mu1 = subs(sol1.mu1, t, delta);
res.mu2 = subs(sol1.mu2, t, delta);
res.delta = solve(res.mu1==res.mu2, delta);
res.delta = eval(res.delta);

求解结果为

与书中结果相同。

于是，得到
$$\{\begin{array}{c}\begin{array}{rl}\lambda (\Delta )& =9\cdot \frac{10}{9}=10,\\ {\mu }_1(\Delta )& =-81\cdot \frac{10}{9}+80=-10,\\ {\mu }_2(\Delta )& =-9\cdot \frac{10}{9}=-10.\end{array}\end{array}\text{\hspace{1em}(12)}$$
在$t=\Delta$处，因$\lambda (\Delta )=10>0$，故
$${\dot{\mu }}_1(\Delta )=9\lambda (\Delta )=90>{\dot{\mu }}_2(\Delta )=\lambda (\Delta )=10.$$
于是在$t=\Delta$的左领域内，有
$$-{\mu }_1(t)>-{\mu }_2(t)\text{\hspace{1em}(13)}$$
由式$(13)$和式$(7)$，有
$${\varphi }_1(t)=1,{\varphi }_2(t)=0,0\le t<\Delta .$$
那么，联立式$(4)$、式$(6)$和式$(8)$，则式$(4)$变为
$$\{\begin{array}{c}\begin{array}{rl}\dot{\lambda }& ={\mu }_1,\\ {\dot{\mu }}_1& =9\lambda ,\\ {\dot{\mu }}_2& =\lambda ,\end{array}\end{array}0\le t\le \Delta$$
考虑终值条件式$(12)$进行求解。

此处编写MATLAB代码求解上述微分方程，仿真代码为

%% 《数理战术学》第7.3节 式（7.33）下面的式子
% 【0 <= t <= delta】 区间
syms delta
eq1 = diff(lambda,t)==mu1;
con1 = lambda(delta)==10;
con2 = mu1(delta)==-10;
con3 = mu2(delta)==-10;
sol2 = dsolve(eq1,eq2,eq3,con1,con2,con3);

求解结果为

书中求解结果为
$$\{\begin{array}{c}\begin{array}{rl}\lambda & =\frac{10}{3}{e}^{3(t-\Delta )}+\frac{20}{3}{e}^{3(\Delta -t)},\\ {\mu }_1& =10e^{3(t-\Delta )}-20e^{3(\Delta -t)},\\ {\mu }_2& =\frac{10}{9}{e}^{3(t-\Delta )}-\frac{20}{9}{e}^{3(\Delta -t)}-\frac{80}{9},\end{array}\end{array}\text{\hspace{1em}(14)}$$
MATLAB求解结果与书中求解结果相同。

通过求解该微分对策问题，发现甲方最优火力分配策略确实发生了改变。在$[0,\Delta )$时间内应集中力量攻击乙方的第一类作战单位，而在$(\Delta ,T]$时间内应集中力量攻击乙方第二类作战单位。

下面来求解相应战斗过程的微分方程组$(1)$的解，并证明甲、乙双方作战单位兵力始终保持如下关系
$$x_1(t)\ge 0,y_1(t)\ge 0,y_2(t)\ge 0$$
下面在MATLAB中编写代码求解微分方程组$(1)$，为方便求解，将式$(1)$再次写在下面
$$\{\begin{array}{c}\begin{array}{rl}{\dot{x}}_1& =-9{\psi }_1{y}_1-{\psi }_2{y}_2,\\ {\dot{y}}_1& =-{\varphi }_1{x}_1,\\ {\dot{y}}_2& =-{\varphi }_2{x}_1.\end{array}\end{array}\text{\hspace{1em}(1)}$$
在$t\in [0,\Delta )$时间内，MATLAB仿真代码为

%% 《数理战术学》第7.3节 求解式（7.19）
% 【0 <= t <= delta】 区间
syms x1(t) y1(t) y2(t)
eq1 = diff(x1)==-9*y1-y2;   con1 = x1(0)==100;
eq2 = diff(y1)==-x1;        con2 = y1(0)==30;
eq3 = diff(y2)==0;          con3 = y2(0)==30;
sol3 = dsolve(eq1,eq2,eq3,con1,con2,con3);

求解结果为

书中求解结果为
$$\{\begin{array}{c}\begin{array}{rl}{x}_1& =100e^{-3t},\\ y_1& =\frac{100}{3}{e}^{-3t}-\frac{10}{3},\\ y_2& =30,\end{array}\end{array}0\le t<\Delta \text{\hspace{1em}(15)}$$
MATLAB求解结果与书中求解结果相同。

在$t\in [\Delta ,T]$时间内，MATLAB仿真代码为

%% 《数理战术学》第7.3节 求解式（7.19）
% 【delta <= t <= T】 区间
delta = 0.3838;
eq1 = diff(x1)==-9*y1-y2;   con1 = x1(0)==subs(sol3.x1,t,delta);
eq2 = diff(y1)==0;          con2 = y1(0)==subs(sol3.y1,t,delta);
eq3 = diff(y2)==-x1;        con3 = y2(0)==subs(sol3.y2,t,delta);
sol4 = dsolve(eq1,eq2,eq3,con1,con2,con3);
syms delta
sol4.x1 = subs(sol4.x1,t,t-delta);
sol4.y1 = subs(sol4.y1,t,t-delta);
sol4.y2 = subs(sol4.y2,t,t-delta);

求解结果为

书中求解结果为
$$\{\begin{array}{c}\begin{array}{rl}{x}_1& =20\sqrt{10}{e}^{\Delta -t}-10\sqrt{10}{e}^{t-\Delta },\\ y_1& =\frac{10}{3}\sqrt{10}-\frac{10}{3},\\ y_2& =20\sqrt{10}{e}^{\Delta -t}+10\sqrt{10}{e}^{t-\Delta }-30\sqrt{10}+30.\end{array}\end{array}\Delta \le t\le T\text{\hspace{1em}(16)}$$
MATLAB求解结果与式$(16)$相比，得到的式子更加复杂。经过核算，

200*exp(-5757/5000) = 63.2388
100*exp(-5757/5000) = 31.6194
300*exp(-5757/5000) = 94.8683
100*exp(-5757/5000)/3-10/3 = 7.2065

$$\{\begin{array}{c}\begin{array}{rl}20\sqrt{10}& =63.2456,\\ 10\sqrt{10}& =31.6228,\\ 30\sqrt{10}& =94.8683,\\ \frac{10}{3}\sqrt{10}-\frac{10}{3}& =7.2076.\end{array}\end{array}$$

两者相差不多，可以认为MATLAB的求解结果与书中结果相同。

至此，使用MATLAB求解微分方程的问题已解决。

将所求得的表达式画图，在MATLAB中编写代码实现，代码为

%% 画图
% 战斗时间
T = (1/6)*log(10) + log(10/9);
delta = (1/6)*log(10);

% 表达式
k = 1;
len = length(0:0.01:T);
x1 = zeros(1,len);
y1 = zeros(1,len);
y2 = zeros(1,len);
for t = 0:0.01:T
    if t <= delta
        x1(k) = 100*exp(-3*t);
        y1(k) = (100/3)*exp(-3*t) - (10/3);
        y2(k) = 30;
        k = k+1;
    else
        x1(k) = 20*sqrt(10)*exp(delta-t) - 10*sqrt(10)*exp(t-delta);
        y1(k) = (10/3)*sqrt(10) - (10/3);
        y2(k) = 20*sqrt(10)*exp(delta-t) + 10*sqrt(10)*exp(t-delta) - 30*sqrt(10) + 30;
        k = k+1;
    end
end
figure('Color',[1,1,1]);box on;
t = 0:0.01:T;
plot(t, x1, 'LineWidth', 1.5);hold on;
plot(t, y1, 'LineWidth', 1.5);
plot(t, y2, 'LineWidth', 1.5);
set(gca,'FontName','Times New Roman',...
        'FontSize',13,...
        'LineWidth',1.3,...
        'XTick',(0:0.05:t(end)),...
        'YTick',(0:10:110))
legend('x1','y1','y2');
xlabel('时间单位','FontName','宋体');ylabel('兵力','FontName','宋体');
axis([0 t(end) 0 110]);

结果为

全部源码如下（作者的MATLAB版本为2018b ）：

clc;clear;
format long;                            % 将数据精度设置为长精度
%% 《数理战术学》第7.3节 式（7.23）、式（7.29）
% 【delta <= t <= T】 区间
syms lambda(t) mu1(t) mu2(t) T
eq1 = diff(lambda,t)==mu2;con1 = lambda(T)==9;
eq2 = diff(mu1,t)==9*lambda;con2 = mu1(T)==-1;
eq3 = diff(mu2,t)==lambda;con3 = mu2(T)==-9;
sol1 = dsolve(eq1,eq2,eq3,con1,con2,con3);
% sol1=dsolve(diff(lambda)==mu2, diff(mu1)==9*lambda, diff(mu2)==lambda,...
%            lambda(T)==9, mu1(T)==-1, mu2(T)==-9);

%% 《数理战术学》第7.3节 式（7.30）
syms delta
T = (1/6)*log(10) + log((10/9));
res.lambda = subs(sol1.lambda, t, delta);
res.mu1 = subs(sol1.mu1, t, delta);
res.mu2 = subs(sol1.mu2, t, delta);
res.delta = solve(res.mu1==res.mu2, delta);
res.delta = eval(res.delta);

%% 《数理战术学》第7.3节 式（7.33）下面的式子
% 【0 <= t <= delta】 区间
syms delta
eq1 = diff(lambda,t)==mu1;
con1 = lambda(delta)==10;
con2 = mu1(delta)==-10;
con3 = mu2(delta)==-10;
sol2 = dsolve(eq1,eq2,eq3,con1,con2,con3);

%% 《数理战术学》第7.3节 求解式（7.19）
% 【0 <= t <= delta】 区间
syms x1(t) y1(t) y2(t)
eq1 = diff(x1)==-9*y1-y2;   con1 = x1(0)==100;
eq2 = diff(y1)==-x1;        con2 = y1(0)==30;
eq3 = diff(y2)==0;          con3 = y2(0)==30;
sol3 = dsolve(eq1,eq2,eq3,con1,con2,con3);

% 【delta <= t <= T】 区间
delta = 0.3838;
eq1 = diff(x1)==-9*y1-y2;   con1 = x1(0)==subs(sol3.x1,t,delta);
eq2 = diff(y1)==0;          con2 = y1(0)==subs(sol3.y1,t,delta);
eq3 = diff(y2)==-x1;        con3 = y2(0)==subs(sol3.y2,t,delta);
sol4 = dsolve(eq1,eq2,eq3,con1,con2,con3);
syms delta
sol4.x1 = subs(sol4.x1,t,t-delta);
sol4.y1 = subs(sol4.y1,t,t-delta);
sol4.y2 = subs(sol4.y2,t,t-delta);

%% 画图
% 战斗时间
T = (1/6)*log(10) + log(10/9);
delta = (1/6)*log(10);

% 表达式
k = 1;
len = length(0:0.01:T);
x1 = zeros(1,len);
y1 = zeros(1,len);
y2 = zeros(1,len);
for t = 0:0.01:T
    if t <= delta
        x1(k) = 100*exp(-3*t);
        y1(k) = (100/3)*exp(-3*t) - (10/3);
        y2(k) = 30;
        k = k+1;
    else
        x1(k) = 20*sqrt(10)*exp(delta-t) - 10*sqrt(10)*exp(t-delta);
        y1(k) = (10/3)*sqrt(10) - (10/3);
        y2(k) = 20*sqrt(10)*exp(delta-t) + 10*sqrt(10)*exp(t-delta) - 30*sqrt(10) + 30;
        k = k+1;
    end
end
figure('Color',[1,1,1]);box on;
t = 0:0.01:T;
plot(t, x1, 'LineWidth', 1.5);hold on;
plot(t, y1, 'LineWidth', 1.5);
plot(t, y2, 'LineWidth', 1.5);
set(gca,'FontName','Times New Roman',...
        'FontSize',13,...
        'LineWidth',1.3,...
        'XTick',(0:0.05:t(end)),...
        'YTick',(0:10:110))
legend('x1','y1','y2');
xlabel('时间单位','FontName','宋体');ylabel('兵力','FontName','宋体');
axis([0 t(end) 0 110]);