# Optimal Control/Runge-Kutta-4

I'm working through Lenhart and Workman's Optimal Control Applied to Biological Models, and I'm trying to apply the Forward-Backward Sweep Method w/ Runge-Kutta-4 as the DE solver to solve

$$\max_{u}\int_{1}^{5}u(t)x(t) - u(t)^{2} - x(t)^{2}\,dt\\\text{subject to}\quad x'(t) = x(t)+u(t),\,x(1) = 2,$$

where both $x$ and $u$ are continuous.

The issue I'm having is that when I solve the problem analytically I get one solution, and when I solve the problem numerically I get an obviously different solution, and I don't know enough about Runge-Kutta Methods to know if I'm programming it incorrectly or if it can even be solved by this method. I'm (fairly) confident in my analytical solution.

I'll first summarize the solution I got analytically.

We form the Hamiltonian $$H = ux - u^{2} - x^{2} + \lambda x + \lambda u,$$ where $\lambda$ is the adjoint. The optimality condition, $$0 = \frac{\partial H}{\partial u} = x - 2u + \lambda \text{ at }u^{*}$$ gives the characterization of the optimal control $$u^{*} = \frac{1}{2}(x+\lambda).$$ The adjoint system is $$\begin{cases} \lambda'(t) = -\frac{\partial H}{\partial x} = -u + 2x - \lambda\\ \lambda(5) = 0\end{cases}.$$ Using the two ODEs and the characterization, we can write down and solve the linear system $$\begin{pmatrix} x\\\lambda \end{pmatrix}' = \begin{pmatrix} 3/2 & 1/2\\3/2 & -3/2 \end{pmatrix}\begin{pmatrix}x\\\lambda \end{pmatrix}.$$ The solution I obtained from this for the state, $x$, is $$x(t) = \frac{6e^{-\sqrt{3}}}{3+2\sqrt{3}+(2\sqrt{3}-3)e^{8\sqrt{3}}}\left(1+\frac{2\sqrt{3}}{3}\right)e^{\sqrt{3}t}-\frac{6e^{9\sqrt{3}}}{3+2\sqrt{3}+(2\sqrt{3}-3)e^{8\sqrt{3}}}\left(1-\frac{2\sqrt{3}}{3}\right)e^{-\sqrt{3}t}.$$ The solutions for the other functions can also be obtained but I'll omit those.

The MATLAB script I wrote to try to solve the problem numerically with the Forward-Backward Sweep Method w/ RK4 is given below. I essentially followed what Lenhart and Workman did for (what seems to me) a similar problem. Note that the requirement for convergence is that $$\delta \|z\| - \|z-oldz\| \geq 0$$ for $z = u,x,\lambda.$ I've only minimally commented my code, if it is not clear I can add more.

test = -1;
delta = 0.01;
N = 4000;
t = linspace(1,5,N+1);
h = 1/N;
h2 = h/2;

u = zeros(1,N+1);
x = zeros(1,N+1);
x(1) = 2;
lambda = zeros(1,N+1);
while(test < 0)
oldu = u;
oldx = x;
oldlambda = lambda;

%RK4 Forward For x
for i = 1:N
k1 = x(i) + u(i);
k2 = x(i) + h2*k1 + 0.5*(u(i) + u(i+1));
k3 = x(i) + h2*k2 + 0.5*(u(i) + u(i+1));
k4 = x(i) + h*k3 + u(i+1);
x(i+1) = x(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end

%RK4 Backwards For lambda
for i = 1:N
j = N + 2 - i
k1 = -u(j) + 2*x(j) - lambda(j);
k2 = -0.5*(u(j) + u(j-1)) - (lambda(j) - h2*k1) + (x(j) + x(j-1));
k3 = -0.5*(u(j) + u(j-1)) - (lambda(j) - h2*k2) + (x(j) + x(j-1));
k4 = -u(j-1) - (lambda(j) - h*k3) + 2*x(j-1);
lambda(j-1) = lambda(j) - (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end

%Update u
u1 = (1/2)*(x + lambda);
u = 0.5*(u1 + oldu);

%Convergence Test
temp1 = delta*sum(abs(u)) - sum(abs(oldu - u));
temp2 = delta*sum(abs(x)) - sum(abs(oldx - x));
temp3 = delta*sum(abs(lambda)) - sum(abs(oldlambda - lambda));
test = min(temp1, min(temp2,temp3));
end

• Other stack overflow project may be able to advise for you. – Takahiro Waki Jan 12 '17 at 3:42

Matlab has a built in implementation of runga-kutta 4. The syntax is as follows:

[timeVar,spaceVar] = ODE45(vectorField,[t_0,t_f],initialCondition)

You should use the built in method to check your implementation for errors. Even if it is correct, it is probably still better to use the included implementation for speed alone.

You can find out more advanced features by typing: help ODE45

• The Dormand-Prince embedded order 4(5) method with adaptive step size is decidedly not the same as the classical RK4 method with fixed step size. Just consult it's Butcher tableau. – Lutz Lehmann Dec 20 '19 at 0:01
• Also, your call variant returning the internal steps with variable step size filled with some interpolation points is unsuitable for the stated purpose of the question. The interpolation provided by the "dense output" variant automatically satisfies the demands of the alternation of forward and backward sweeps where one function is input for the integration of the other. – Lutz Lehmann Dec 20 '19 at 0:11