Example 3¶
Corresponds to Example 5.1 in [1].
[1]:
import pymloc
import numpy as np
import jax.numpy as jnp
import warnings
warnings.filterwarnings('ignore')
Creating the variables object¶
The variables for the different levels are defined as follows.
[2]:
from pymloc.model.variables import InputStateVariables
from pymloc.model.variables import NullVariables
from pymloc.model.variables import ParameterContainer
from pymloc.model.variables.time_function import Time
from pymloc.model.domains import RNDomain
loc_vars = InputStateVariables(1, 1, time=Time(0., 2.))
hl_vars = ParameterContainer(1, domain=RNDomain(1))
variables2 = (hl_vars, loc_vars)
ll_vars = NullVariables()
Creating the parameter optimal control object¶
First, we need to create a parameter dependent optimal control problem.
We need objective and constraint.
Creating the control system¶
The parameter dependent control system is defined by
[3]:
from pymloc.model.control_system.parameter_dae import LinearParameterControlSystem
def e(p, t):
return jnp.array([[1.]])
def a(p, t):
return jnp.array([[-1.]])
def b(p, t):
return jnp.array([[1.]])
def c(p, t):
return jnp.array([[1.]])
def d(p, t):
return jnp.array([[0.]])
def f(p, t):
return jnp.array([0.])
param_control = LinearParameterControlSystem(ll_vars, *variables2, e, a, b, c,
d, f)
Creating the constraint object¶
[4]:
from pymloc.model.optimization.parameter_optimal_control import ParameterLQRConstraint
def initial_value(p):
return jnp.array([2.])
time = Time(0., 2.)
pdoc_constraint = ParameterLQRConstraint(*variables2, param_control,
initial_value)
Creating the objective function¶
The objective function is defined by
[5]:
from pymloc.model.optimization.parameter_optimal_control import ParameterLQRObjective
def q(p, t):
return jnp.array([[p**2. - 1.]])
def s(p, t):
return jnp.zeros((1, 1))
def r(p, t):
return jnp.array([[1.]])
def m(p):
return jnp.array([[0.]])
time = Time(0., 2.)
pdoc_objective = ParameterLQRObjective(*variables2, time, q, s, r, m)
Create the parameter dependent optimal control object¶
[6]:
from pymloc.model.optimization.parameter_optimal_control import ParameterDependentOptimalControl
pdoc_object = ParameterDependentOptimalControl(*variables2, pdoc_objective,
pdoc_constraint)
The neccessary conditions can be obtained as follows
[7]:
neccessary_conditions = pdoc_object.get_bvp()
e = neccessary_conditions.dynamical_system.e(2., 3.)
a = neccessary_conditions.dynamical_system.a(2., 3.)
print("E =\n {},\nA =\n {}".format(e, a))
E =
[[ 0. 1. 0.]
[-1. 0. 0.]
[-0. 0. 0.]],
A =
[[ 0. -1. 1.]
[-1. 3. 0.]
[ 1. 0. 1.]]
Setting up the nonlinear least squares problem¶
Define the residual function¶
We define the residual function by using the true solution for the parameter value at 2 and comparing it with the current solution for the current parameter value.
[8]:
def compute_ref_sol(theta, time, x0, t):
tf = time.t_f
t0 = time.t_0
exp0 = np.exp(2 * theta * (tf - t0))
exp1 = np.exp(-(t + t0) * theta)
exp2 = np.exp(2 * t * theta)
exp3 = np.exp(2 * tf * theta)
tmp1 = theta + exp0 * (theta + 1) - 1
tmp2 = np.array([
-(exp2 - exp3) * (theta**2 - 1),
(exp2 * (theta - 1) + exp3 * (theta + 1)),
(exp2 - exp3) * (theta**2 - 1)
])
refsol = tmp1**-1 * tmp2 * exp1 * x0
return refsol
[9]:
theta = 2.
t0 = 0.
tf = 2.
x01 = 2.
time = Time(t0, tf)
t2 = 1.3
t1 = 1.
def f_nlsq(ll_vars, hl_vars, loc_vars):
sol1 = compute_ref_sol(theta, time, x01, t1)
sol2 = compute_ref_sol(theta, time, x01, t2)
sol0 = compute_ref_sol(theta, time, x01, t0)
solf = compute_ref_sol(theta, time, x01, tf)
f1 = ll_vars(t1)[1] - sol1[1]
f2 = ll_vars(t2)[1] - sol2[1]
f0 = ll_vars(t0)[1] - sol0[1]
ff = ll_vars(tf)[1] - solf[1]
return np.hstack((f0, f1, f2, ff))
Define the NonlinearLeastSquares optimization object¶
[10]:
import pymloc.model.optimization as optimization
variables = (variables2[1], NullVariables(), variables2[0])
nlsq_obj = optimization.objectives.NonLinearLeastSquares(*variables, f_nlsq)
nlsq = optimization.NonLinearLeastSquares(nlsq_obj, *variables)
Setting up the bilevel optimal control problem¶
Compute solution¶
We are now able to set up the multilevel optimal control problem by setting the optimizations and variables in the corresponding order
[11]:
import logging
logger = logging.getLogger(__name__)
logging.getLogger().setLevel(logging.INFO)
optimizations = [nlsq, pdoc_object]
variables = (variables2[0], variables2[1])
variables[0].current_values = np.array([1.])
variables[1].current_values = np.array([])
variables[1].time.grid = np.array([1., 1.3])
We also only want sensitivities of the x component and thus, set the selector accordingly.
[12]:
pdoc_object.ll_sens_selector_shape = (1, 3)
pdoc_object.ll_sens_selector = lambda p: np.array([[0., 1., 0.]])
Then, wen can initialize and run the multilevel optimal control problem
[13]:
from pymloc import MultiLevelOptimalControl
logger = logging.getLogger("pymloc.solvers.nonlinear.gauss_newton")
logger.setLevel(logging.DEBUG)
logging.getLogger().handlers[0].filters[0].__class__.max_level = 3
mloc = MultiLevelOptimalControl(optimizations, variables)
np.set_printoptions(precision=8)
mloc.init_solver(abs_tol=1e-6, rel_tol=1e-6)
mloc.highest_opt.local_level_variables.associated_problem.solver_instance.upper_eta = 0.1
solution = mloc.solve()
logger.info("Solution: {}".format(solution.solution))
Starting solver MultiLevelIterativeSolver
Current option values:
abs_tol: 1e-06
rel_tol: 1e-06
max_iter: 10
Starting solver GaussNewton
Current option values:
abs_tol: 1e-06
rel_tol: 1e-06
max_iter: 20
Starting iteration: 0
Updating lower level variables...
Current residual: [-0.66603219]
Current allowed lower level tolerance: 0.023322921822300516
New x value:
[1.62565779]
New jac value:
[[ 0. ]
[-0.69271902]
[-0.64637977]
[-0.40849072]]
Starting iteration: 1
Updating lower level variables...
Current residual: [-0.09531969]
Current allowed lower level tolerance: 0.008435398848567852
New x value:
[1.92942009]
New jac value:
[[ 0. ]
[-0.4116986 ]
[-0.33576026]
[-0.17766904]]
Starting iteration: 2
Updating lower level variables...
Current residual: [-0.01137531]
Current allowed lower level tolerance: 0.0019700405466293908
New x value:
[2.00504389]
New jac value:
[[ 0. ]
[-0.29951971]
[-0.22281496]
[-0.10517187]]
Starting iteration: 3
Updating lower level variables...
Current residual: [0.00027421]
Current allowed lower level tolerance: 5.913012890989498e-05
New x value:
[2.00281763]
New jac value:
[[ 0. ]
[-0.27465787]
[-0.19906282]
[-0.09003094]]
Starting iteration: 4
Updating lower level variables...
Current residual: [0.00033822]
Current allowed lower level tolerance: 7.339535863148346e-05
New x value:
[2.00005288]
New jac value:
[[ 0. ]
[-0.27377794]
[-0.19853268]
[-0.08924611]]
Starting iteration: 5
Updating lower level variables...
Current residual: [-1.39727571e-06]
Current allowed lower level tolerance: 3.036733193512886e-07
New x value:
[2.00006423]
New jac value:
[[ 0. ]
[-0.27453965]
[-0.19927597]
[-0.08970326]]
Starting iteration: 6
Updating lower level variables...
Current residual: [7.91042637e-06]
Current allowed lower level tolerance: 1.7204548910448125e-06
New x value:
[1.99999989]
New jac value:
[[ 0. ]
[-0.27437847]
[-0.19913989]
[-0.08950784]]
Starting iteration: 7
Updating lower level variables...
Current residual: [-3.03250968e-07]
Current allowed lower level tolerance: 6.595525325505868e-08
Solution: [2.00000235]