# Problem with solving an optimal control system with Shooting Method

I try to solve a system of differential equation by shooting method. I looked at questions asked on the site but I could not find that could help me.

My functions are ;

u[c_] := (c^(1 - σ) - 1^(1 - σ))/(1 - σ)
f[s_] := s + g s (1 - s/sbar1)
h[s_] := (2 hbar)/(1 + Exp[η (s/sbar - 1)])
co[a_] := ϕ  (a^2)/2
ψ[k_] := wbar (ω + (1 - ω) Exp[-γ k])


Calibration I use

paramFinal = {σ -> 1.5, ρ -> 0.025, g -> 0.05, sbar -> 51.25, η -> 3.5, hbar -> 0.5,
χ -> 10, ω -> 0.1, δ -> 0.065, ϕ -> 0.05, sbar1 -> 5, wbar -> 1, γ -> 0.6};


Differential system that I try to solve is

dec1 = c'[t] == -(u'[c[t]]/u''[c[t]]) (f'[s[t]] - (ρ + h[s[t]]) -
h'[s[t]]/(ρ + h[s[t]]) (1/u'[c[t]] (u[c[t]] - ψ[k[t]] h[s[t]] - co[a[t]]))
- (ψ[k[t]]  h'[s[t]])/u'[c[t]] - h'[s[t]]/(ρ + h[s[t]]) (f[s[t]] - c[t])
- (h'[s[t]] co'[a[t]] )/((ρ + h[s[t]]) u'[c[t]]) (a[t] - δ k[t]));

dea1 = a'[t] == (co'[a[t]]/co''[a[t]] (ρ + h[s[t]] + (ψ'[k[t]] h[s[t]])/co'[a[t]] + δ));

des1 = s'[t] == (f[s[t]] - c[t]);

dek1 = k'[t] == (a[t] - δ k[t]);


Variables $a$ and $c$ are control, $s$ and $k$ are state variables. I implement the following code for Shooting Method. So, I start with a guess, in order to find optimal initial values for control variables.

ODEs2[ω1_, ω2_] := {dec1, dea1, des1, dek1, c[0] == ω1, a[0] == ω2, k[0] == 4.1,
s[0] == 95}


Initial conditions for state variables are given.

Soln2[ω1_, ω2_] := NDSolve[ODEs2[ω1, ω2] /. paramFinal, {c, a, s, k}, {t, 0, 0.5}]


I put an end condition ;

EndCondition[ω1_?NumericQ, ω2_?NumericQ] := {First[c[t] /. Soln2[ω1, ω2]] /.
t -> 0.5, First[a[t] /. Soln2[ω1, ω2]] /. t -> 0.5}


with some arbitrary values ;

EndCondition[7, 0.2]


It is sure that this first guess is not correct.

By FindRoot, I write ;

initVals2 = FindRoot[EndCondition[α1, β1] == 0, {α1, .01, 100}, {β1, 0.01, 200}]


The problem is that I always get a stiffness problem as follows ;

NDSolve::ndsz: At t == 0.012993888073465564, step size is effectively zero; singularity or stiff system suspected. >>

I also get other similar error messages.

What can I do in order to solve this issue ? I looked at some options like EventLocator in order to stop integration at $t$ where there is stiffness problem but it did not help. I really appreciate any help and hints.

In general, when FindRoot is initialized with two initial guesses for each variable, it uses them to estimate how the function to be solved varies locally with those variables. Giving initial guesses that are far apart, as in the question, defeats this purpose. Trying something like

initVals2 = FindRoot[EndCondition[α1, β1] == 0, {α1, 6, 8}, {β1, 0.01, 0.03}] // Chop
(* {α1 -> 0, β1 -> 0.0194729} *)


avoids the error cited in the question.

Simpler Solution

Define

ode = {Expand[dec1], dea1, des1, dek1} /. paramFinal;


(Expand is used to eliminate a spurious singularity at c = 0.)

sol = NDSolve[{ode, c[1/2] == 0, a[1/2] == 0, k[0] == 41/10, s[0] == 95},
{c, a, s, k}, {t, 0, 1/2}];


from which the initial conditions on c and a can be obtained from

{(c[0] /. sol), a[0] /. sol} // Chop
(* {{-2.66636*10^-7}, {0.0194729}} *)


In fact, c actually is equal to zero, as can be seen by evaluating

ode[[1]] /. c[t] -> 0
(* Derivative[1][c][t] == 0 *)


If desired, the remaining equations can be solved by

sol1 = NDSolve[{ode[[2 ;; 4]] /. c[t] -> 0, a[1/2] == 0, k[0] == 41/10, s[0] == 95},
{a, s, k}, {t, 0, 1/2}];
a[0] /. sol1
(* {0.0194729} *)

• Interesting…I'd like to mention that, v10 seems to be necessary to solve the equation with the simpler solution, in v9.0.1 and v8.0.4 it failed. And the Expand seems to be redundant, at least according to my test on Wolfram Cloud. BTW, how did you find the good initial guess for the not-that-simple solution? – xzczd Feb 14 '16 at 3:24
• @xzczd I had intended to use a guess based on the line of code in the question, EndCondition[7, 0.2] but mistyped the second parameter. With respect to Expand, I found that dec1 /. c[t] -> 0` generated an error otherwise. – bbgodfrey Feb 14 '16 at 4:57
• @bbgodfrey I also could not run the simpler solution on v9.0.1, thanks so much for your help. – optimal control Feb 14 '16 at 14:33