本文主要基于以下参考：

[1] John T. Betts. Survey of Numerical Methods for Trajectory Optimization.

[2] Anil V. Rao. A Survey of Numerical Methods For Optimal Control.

[3] John T. Betts. Practical Methods for Optimal Control and Estimation Using Nonlinear Programming 2nd.

[4] A E. Bryson. Applied Optimal Control.

[5] KIRK. Optimal Control Theory: An Introduction.

[6] Matthew Kelly. An Introduction to Trajectory Optimization: How to Do Your Own Direct Collocation.

本节不涉及具体的理论知识，主要对前面的部分内容进行回顾整理。

在第9节中介绍了 Zermelo 问题，并对其最优必要条件进行了推导，不过那里的公式大部分都是直接从书里直接摘抄过来的，并没有重新推导，本节便对这些公式进行推导，并在最后通过GPOPS对该问题进行求解。

第一部分

概述：主要对第9节中的控制律  的表达式进行了推导。

首先对问题重新叙述(Zermelo’s problem)：考虑一艘船必须穿过一片有较强水流的区域，该水流的幅值和方向为位置的已知函数：

其中，  为平面坐标，  为水流在  和  方向的速度分量，船体相对于水流的速度为  (常数)。该最优控制问题是指控制该船使得在最短时间内从  点到达  点。

运动方程如下：

其中，  为船头与坐标系之间的夹角，  代表船在水中的位置。

系统哈密顿函数为：

于是，最优必要条件如下：

此外，由于哈密顿函数(10)不是时间的显式函数，因此哈密顿在整个过程中为常数。此外，根据之前的最优条件：

因此，哈密顿函数必须为零，即  ，因此根据式(3)式(4)后一式可以解得：

注意，上述公式的推导仅使用了两个必要条件，此时将得到的协态代入式(4)的前两式，可以得到：

其中，  ，对上式进行整理，将带  的整理到一起，另外的在一起，可以得到：

继续对左边进行整理，将  带入以及  带入整理可得：

于是可以得到：

第二部分

概述：主要对第10节中的最后的公式(26)进行推导。

问题描述参考第10节，列出下述公式：

显然，此时在这两个区域中，  和  都是常数，但是在条件

处不连续跳变，由公式：

可以得到：

同样应用公式

到本问题中可以得到：

由最优控制条件可得：

由轨迹的连续性，我们有  以及  。因此，最优路径由两条直线构成，因此可得：

以上这些方程构成了含有8个未知数的方程组，下面考虑求解这些方程，因为我们只对  和

 感兴趣，因此首先将(16)代入到式(15)中便可以得到：

由条件(14)可得，并约去  ：

然后进行整理就可以得到：

注意，上式是本节出现的原因，我发现书上给的公式符号有问题，书上是负号，但实际上应该是正号，这一点已经通过数值进行了验证。

然后同样，可以首先由式(17)解出  ，然后代入到第二式子中，便可以得到：

于是，在求得这两个解以后，可以依次计算得到

第三部分数值仿真

例子1

首先对第9节的问题进行仿真，得到的仿真结果如下：

显然，数值仿真结果与之前采用间接法计算得到的解的形态是相同的，然而直接法的初值收敛范围很大，因此即使初值较差，也可以收敛到相应的最优解。

例子2

第10节的右端函数不连续情形

显然，仿真结果也验证了前面解析推导的结果，最优轨迹由两段直线构成，以  为分界线，最优控制律在两段分别为常数。

进一步我们采用前面推导得到的方程，并结合fsolve来对最优解进行求解，计算代码如下：

% fsolve

x0  = [0.7962; 1.0260];
[x, feval] = fsolve(@fun,x0);

%%
xneg = x(1);    xpos = x(2);
V = 2;

lambdax = -1/(V*cos(xneg)+tan(xneg)*V*sin(xneg));

lambdayneg = -(1+lambdax*V*cos(xneg))/V/sin(xneg);

lambdaypos = -(1+lambdax*(V*cos(xpos)+0.6*V))/V/sin(xpos);

%%

function ret=fun(x)
    epsilon = 0.6;
    a = 1;
    b = 1;
    xneg = x(1);    xpos = x(2);
    ret(1) = sec(xneg)-sec(xpos)-epsilon;
    ret(2) = cot(xneg)-a+b*(cot(xpos)+epsilon*csc(xpos));
end

计算结果如下:
x =

    1.3905
    1.3684

与数值计算得到的结果相同，下面进一步对系统协态变量进行计算，由已知系统哈密顿  ，因此由：

此时，  已知，因此上式中一共有三个未知数四个方程，可以计算得到计算公式：

计算得到：

然后由下式可以计算出  ：

以上的计算结果与数值计算结果相吻合。

总结

本节主要对前面最短时间航行问题中的几个公式进行梳理，并对关键公式进行了推导，发现了书中公式的符号存在问题。然后通过数值仿真实验与解析得到的计算公式进行了验证。