Problem with solving an optimal control system with Shooting Method

Question

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.

Answer 1

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} *)