I am solving a sensitivity analysis problem, the direct equation corresponds to an initial-value problem, and the adjoint problem is a boundary value problem. The initial value problem is solved (perfectly) with scipy.integrate.odeint() and the boundary value problem should be solved with scipy.integrate.solve_bvp(). I will share the 4 routines I did, to solve these problems and point later the execution error in a particular example.
ODEINT STUFF
def ode_sys(sol,t,A,B,D):
"""
This function returns the order reduction for the ODE's (1) and (2)
from Singh2020.pdf
A , B, D: are (known) scalar parameters
x_sol , y_sol: are the vocal fold displacements
u_sol , v_sol: are the vocal fold velocities
You need to convert the 2-second order ODEs into 4-first order ODEs
"""
x_sol , u_sol , y_sol , v_sol = sol
x_der = u_sol
u_der = (A - B * (1 + np.power(x_sol,2))) * u_sol + (0.5 *D-1) * x_sol + A * v_sol
y_der = v_sol
v_der = (A - B * (1 + np.power(y_sol,2))) * v_sol - (0.5 *D+1) * y_sol + A * u_sol
return [x_der,u_der,y_der,v_der]
def ode_solver(A,B,D,sol_0,t):
"""
This function integrates the coupled ODEs that represent the Van der Pol oscillators
A :: alpha
B :: beta
D :: delta
sol_0 :: initial condition [x(0),x'(0),y(0),y'(0)]
t :: time array
u0 :: computed glottal flow
sol :: displacements of both vocal folds vs time
"""
# Here I define constants of the model (physiological props.)
c_til,d_len,x_0 = physical_props() # they are scalar values
sol = odeint(ode_sys,sol_0,t,args=(A,B,D))
u0 = c_til * d_len * (2 * x_0 + sol[:,0] + sol[:,2])
return u0,sol
SOLVE_BVP STUFF
def adj_sys(adj_sol,t,voc_sol,data,A,B,D):
"""
This function returns the order reduction for ODE's 20-23 from Singh2020B.pdf
A , B, D: are (known) scalar parameters
x_adj , y_adj: are the adjoint displacements
u_adj , v_adj: are the adjoint velocities
x_sol , y_sol: are the vocal fold displacements (known data)
u_sol , v_sol: are the vocal fold velocities (known data)
You need to convert the 2-second order ODEs into 4-first order ODEs
"""
print('inside sys,adj',adj_sol.shape)
print('inside sys,voc',voc_sol.shape)
x_adj , u_adj , y_adj , v_adj = adj_sol
x_sol , u_sol , y_sol , v_sol = voc_sol
# Here I define constants of the model (physiological props.)
c_til,d_len,x_0 = physical_props()
u0 = c_til * d_len * (2 * x_0 + x_sol + y_sol)
R = u0 - data
x_adj_der = u_adj
u_adj_der = A * (x_adj + y_adj) - 2.0*c_til*d_len*R - B * (1.0 + np.power(x_sol,2)) * x_adj + (0.5*D - 1.0 - 2.0*B*x_sol*u_sol)
y_adj_der = v_adj
v_adj_der = A * (x_adj + y_adj) - 2.0*c_til*d_len*R - B * (1.0 + np.power(y_sol,2)) * y_adj - (0.5*D + 1.0 + 2.0*B*y_sol*v_sol)
return [x_adj_der,u_adj_der,y_adj_der,v_adj_der]
def bc(x_adj , u_adj , y_adj , v_adj):
"""
This function returns the boundary conditions from Eqs.24-27 from Singh2020B.pdf
"""
# These return values are what we want to be 0:
return x_adj[-1] , u_adj[-1] , y_adj[-1] , v_adj[-1]
def adj_solver(A,B,D,voc_sol,t):
"""
This function integrates the coupled ODEs that represent the adjoint state for Van der Pol oscillators
A :: alpha
B :: beta
D :: delta
voc_sol :: Vocal fold displacements and velocities
t :: time array
"""
#result = solve_bvp(adj_sys, bc, t, voc_sol, voc_sol , data , A , B , D)
result = solve_bvp(lambda y , t: adj_sys(y,t,voc_sol,data,A,B,D),
lambda y0, y1, y2, y3: bc(y0, y1, y2, y3),
t, voc_sol)
adj_sol = [result.y[0],result.y[1],result.y[2],result.y[3]]
return adj_sol
As you can see to solve the adjoint problem I first need to solve the initial value problem. And the solutions of the initial value problem are inputs for the adjoint problem...
The content of the main function
# Here I define constants of the model (physiological props.)
c_til,d_len,x_0 = physical_props()
sol_0 = [0,0.1,0,0.1] # Initial state from papers (reference) (x(0),x'(0),y(0),y'(0))
theta = [c_til,d_len,x_0]
# Here I define a synthetic dataset just to test the integration
D = 0.3
A = 0.6 * D
B = 0.2
t = np.linspace(0,60*pi,1000)
data,voc_sol = ode_solver(A,B,D,sol_0,t)
# Solve the adjoint state
print(voc_sol.T.shape)
adj_state = adj_solver(A,B,D,voc_sol.T,t)
So to summarize ode_solver() works perfectly, the results are as expected. But adj_solver() fails to produce the solution showing the following error message:
(4, 1000)
inside sys,adj (1000,)
inside sys,voc (4, 1000)
Traceback (most recent call last):
File "demo.py", line 166, in <module>
adj_state = adj_solver(A,B,D,voc_sol.T,t)
File "demo.py", line 123, in adj_solver
t, voc_sol)
File "/home/valdez/anaconda3/lib/python3.7/site-packages/scipy/integrate/_bvp.py", line 1063, in solve_bvp
f = fun_wrapped(x, y, p)
File "/home/valdez/anaconda3/lib/python3.7/site-packages/scipy/integrate/_bvp.py", line 648, in fun_p
return np.asarray(fun(x, y), dtype)
File "demo.py", line 121, in <lambda>
result = solve_bvp(lambda y , t: adj_sys(y,t,voc_sol,data,A,B,D),
File "demo.py", line 87, in adj_sys
x_adj , u_adj , y_adj , v_adj = adj_sol
ValueError: too many values to unpack (expected 4)
See in the terminal prints, the arrays for the adjoint solver are wrong... Even if I supplied the initial guess with the correct shape array.
Any help, advice, and suggestion will be very well-acknowledged, will save my day.
I forgot to add the answers I used to write the adjoint solver Answer 1 and Answer 2 Both of them used solve_bvp()... In my case I am having trouble making a callable function without extra arguments.
