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

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

在《数理战术学》第79页提到：

最优火力分配策略随时间变化的根本原因在于相对作战指数向量与实际交战过程中各类作战单位的作战效果很不一致，改变例1中乙方作战指数向量，则最优火力分配策略会发生相应变化，请读者自行证明。

此文为《数理战术学》该思考题的证明。

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

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

$$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}9\\ 1\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)-9y_1(T)-y_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=-9,\\ {\mu }_2(T)& =\frac{\partial \Theta }{\partial y_2}{|}_T=-1\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)$时（即变量${\mu }_1(t)$让$H$负得更少），有
$${\varphi }_1(t)=1,{\varphi }_2(t)=0.\text{\hspace{1em}(7)}$$
当${\mu }_1(t)>{\mu }_2(t)$时（即变量${\mu }_2(t)$让$H$负得更少），有
$${\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)}$$
由式$(7)$有
$${\varphi }_1(t)=1,{\varphi }_2(t)=0,\Delta <t\le T$$
那么，联立式$(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}\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);

因在$t=T$处有${\mu }_1(t)<{\mu }_2(t)$，此处需要考察${\mu }_1(t)$和${\mu }_2(t)$的初值。

如果有${\mu }_1(0)<{\mu }_2(0)$，则在$t\in [0,T]$内恒有${\mu }_1(t)<{\mu }_2(t)$。

如果有${\mu }_1(0)>{\mu }_2(0)$，则在$t\in [0,\Delta ]$内恒有${\mu }_1(t)<{\mu }_2(t)$，在$t\in (\Delta ,T]$内恒有${\mu }_1(t)>{\mu }_2(t)$。

即如下图（A）、（B）所示情况。

上述MATLAB代码已经求得${\mu }_1(t)$和${\mu }_2(t)$的表达式，因此只需再在MATLAB中编写代码求解${\mu }_1(0)$和${\mu }_2(0)$即可。

MATLAB代码为

%% 《数理战术学》第7.3节 式（7.30）
syms delta
res.lambda = subs(sol1.lambda, t, 0);
res.mu1 = subs(sol1.mu1, t, 0);
res.mu2 = subs(sol1.mu2, t, 0);
if eval(res.mu1) < eval(res.mu2)
   disp('mu1 恒小于 mu2'); 
else
    disp('mu1 不恒小于 mu2'); 
end

结果为

eval(res.mu1)=-52.101934085986272
eval(res.mu2)=-14.677992676220697

说明${\mu }_1(0)<{\mu }_2(0)$，那么有$t\in [0,T]$内${\mu }_1(t)<{\mu }_2(t)$。

即如下图所示情况。

也就是说，在$t\in [0,T]$内，甲方的火力分配策略不会改变，持续攻击乙方第一类作战单位。该结果与上一篇文章的结果有所区别。

证明了书中结论：作战指数改变后，会导致最优火力分配策略改变。

此处可以看出，甲方与乙方第一类作战单位的交战强度比甲方与乙方第二类作战单位的交战强度更加“激烈”，并且考虑乙方第一类作战单位的作战指数大于乙方第二类作战单位的作战指数，说明乙方第一类作战单位的重要程度强于第二类作战单位的重要程度。

按照数理战术学理论，双方倾向于先完成交战强度“激烈”的交战格局。因此甲方先与乙方第一类作战单位是合理的。在交战过程中，先消灭交战强度“激烈”的敌方作战单位，减少对己方的威胁，保存己方实力，再消灭对己方威胁程度较低的敌方作战单位，符合常规作战认知。