本文主要基于以下参考：

[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.

最优控制问题数值方法-间接法/变分法 (7

理论部分

本部分涉及到的内容为[4]中的3.5节部分，考虑最优控制问题：

此外，考虑内点约束：

类似之前，考虑到关系式：

将式(15)带入到式(14)中，可以得到：

因此，可以得到最优必要条件如下：

至此，我们得到了一个三点边值问题。

下面给出一个具体的例子：

最优控制问题(Interior point constraint problem)：

满足边界约束：

内点约束：

类似前面，首先考虑写出最优必要条件：

鉴于这个问题的系统比较简单，因此可以不借助于编程实现：

首先将最优控制 u 求得：

代入到系统动力学中：

然后可以对上式继续微分，并将协态代入得到：

上式的解如下：

现在我们面对的是一个三点边值问题，通过选择 π 来使得条件满足，首先对于第一段可以得：

根据上式可以计算得到：

进一步根据式(26)，仅将常数改变，得：

由关系式：

可以解得相应得数值，进一步得到协态函数值：

可以计算得到相应的：

因此便可以计算得到相应的数值 π ，即：

下面给出用MATLAB计算的代码示例：

A = [1 1;exp(0.75) exp(-0.75)];
b = [1;0.9];

ret = A \ b;
C1 = ret(1);    C2 = ret(2);

lambda_neg=-2*(C1*exp(0.75)-C2*exp(-0.75));

% then compute C3 C4
A1 = [exp(0.75) exp(-0.75);exp(1) exp(-1)];
b1 = [0.9;0.75];

ret1 = A1 \ b1;
C3 = ret1(1);    C4 = ret1(2);

lambda_pos=-2*(C3*exp(0.75)-C4*exp(-0.75));

lambda_pi = lambda_neg - lambda_pos;

% state
t1 = linspace(0, 0.75);
s1 = C1*exp(t1)+C2*exp(-t1);
t2 = linspace(0.75, 1);
s2 = C3*exp(t2)+C4*exp(-t2);
figure(1); clf
plot(t1,s1,'r-',t2,s2,'g-','LineWidth',2)
legend('state seg1','state seg2','Location','NorthEast')
grid on;
%axis([0 0.9 0 7])
%title('HBVP Solution')
xlabel('Time');ylabel('States')

% costate
t1 = linspace(0, 0.75);
cs1 = -2*(C1*exp(t1)-C2*exp(-t1));
t2 = linspace(0.75, 1);
cs2 = -2*(C3*exp(t2)-C4*exp(-t2));
figure(2); clf
plot(t1,cs1,'r-',t2,cs2,'g-','LineWidth',2)
hold on;
plot([t1(end) t1(end)],[cs1(end) cs2(1)], 'k--','LineWidth',2)
legend('costate seg1','costate seg2','Location','NorthWest')
grid on;
%axis([0 0.9 0 7])
%title('HBVP Solution')
xlabel('Time');ylabel('Costates')

% control
u1 = -0.5*cs1;
u2 = -0.5*cs2;
figure(3); clf
plot(t1,u1,'r-',t2,u2,'g-','LineWidth',2)
hold on;
plot([t1(end) t1(end)],[u1(end) u2(1)], 'k--','LineWidth',2)
legend('control seg1','control seg2','Location','NorthWest')
grid on;
%axis([0 0.9 0 7])
%title('HBVP Solution')
xlabel('Time');ylabel('Control')

% performance values
disp('the performance values is: ');
trapz([t1 t2], [s1.^2+u1.^2 s2.^2+u2.^2])

不过，MATLAB中的bvp4c也可以用来求解多点边值问题，因此也可以用来求解该问题，最后从求解得到的协态值计算得到 [公式] 的值。

有兴趣的读者可以尝试这种方法并与解析推导得到的结果进行比对。