Example 3¶
Corresponds to Example 5.1 in [Ban2020_ParameterDependentMultilevel].
[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)
gauss_newton = mloc.highest_opt.local_level_variables.associated_problem.solver_instance
gauss_newton.upper_eta = 0.01
gauss_newton.save_intermediate = True
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 ||J^T r||_2: 0.6660321871331705
Current ||r||_2: 0.6472009202931233
Current ||J||_2: 1.0317611699564475
Current cond(J): 1.0
Current allowed lower level tolerance: 0.0025501543923860965
New x value:
[1.62565779]
New jac value:
[[ 0. ]
[-0.69271902]
[-0.64637977]
[-0.40849072]]
Starting iteration: 1
Updating lower level variables...
Current residual ||J^T r||_2: 0.09021829248969856
Current ||r||_2: 0.16544194785738267
Current ||J||_2: 0.5456313369531584
Current cond(J): 1.0
Current allowed lower level tolerance: 0.0009145542277046758
New x value:
[1.92869503]
New jac value:
[[ 0. ]
[-0.40312896]
[-0.32529563]
[-0.17141571]]
Starting iteration: 2
Updating lower level variables...
Current residual ||J^T r||_2: 0.010023113613176841
Current ||r||_2: 0.026234802667399967
Current ||J||_2: 0.38206638166805335
Current cond(J): 1.0
Current allowed lower level tolerance: 0.00020425261292175866
New x value:
[1.99735838]
New jac value:
[[ 0. ]
[-0.29575558]
[-0.21925259]
[-0.10213549]]
Starting iteration: 3
Updating lower level variables...
Current residual ||J^T r||_2: 0.00033601046465835875
Current ||r||_2: 0.000953972969659078
Current ||J||_2: 0.35227148193310337
Current cond(J): 1.0
Current allowed lower level tolerance: 8.387938329980655e-06
New x value:
[2.00006606]
New jac value:
[[ 0. ]
[-0.27553411]
[-0.20008763]
[-0.090228 ]]
Starting iteration: 4
Updating lower level variables...
Current residual ||J^T r||_2: 7.553927904696164e-06
Current ||r||_2: 2.15665120392612e-05
Current ||J||_2: 0.35068348973442637
Current cond(J): 1.0
Current allowed lower level tolerance: 1.9041429186991604e-07
New x value:
[2.00000464]
New jac value:
[[ 0. ]
[-0.27440401]
[-0.19915694]
[-0.0895425 ]]
Starting iteration: 5
Updating lower level variables...
Current residual ||J^T r||_2: 5.705364273345304e-07
Current ||r||_2: 1.6437829701113807e-06
Current ||J||_2: 0.3506677917787622
Current cond(J): 1.0
Current allowed lower level tolerance: 1.4383951618479082e-08
Solution: [2.]
Create table of intermediate data
[14]:
import tabulate
saved = solution.params['intermediates']
print(tabulate.tabulate(saved,headers='keys',tablefmt='latex_booktabs',floatfmt='.7e'))
\begin{tabular}{rrrrrrr}
\toprule
iteration & x & r & J & JTr & atol & cond \\
\midrule
0 & 1.0000000e+00 & 6.4720092e-01 & 1.0317612e+00 & 6.6603219e-01 & 1.0000000e-06 & 1.0000000e+00 \\
1 & 1.6256578e+00 & 1.6544195e-01 & 5.4563134e-01 & 9.0218292e-02 & 2.5501544e-03 & 1.0000000e+00 \\
2 & 1.9286950e+00 & 2.6234803e-02 & 3.8206638e-01 & 1.0023114e-02 & 9.1455423e-04 & 1.0000000e+00 \\
3 & 1.9973584e+00 & 9.5397297e-04 & 3.5227148e-01 & 3.3601046e-04 & 2.0425261e-04 & 1.0000000e+00 \\
4 & 2.0000661e+00 & 2.1566512e-05 & 3.5068349e-01 & 7.5539279e-06 & 8.3879383e-06 & 1.0000000e+00 \\
5 & 2.0000046e+00 & 1.6437830e-06 & 3.5066779e-01 & 5.7053643e-07 & 1.9041429e-07 & 1.0000000e+00 \\
& 2.0000000e+00 & & & & 1.4383952e-08 & \\
\bottomrule
\end{tabular}
Perform the same computation for lower tolerances.¶
[15]:
variables[0].current_values = np.array([1.])
variables =(variables[0], InputStateVariables(1, 1, time=Time(0., 2.)))
variables[1].current_values = np.array([])
variables[1].time.grid = np.array([1., 1.3])
mloc = MultiLevelOptimalControl(optimizations, variables)
mloc.init_solver(abs_tol=1e-1, rel_tol=1e-1)
gauss_newton = mloc.highest_opt.local_level_variables.associated_problem.solver_instance
gauss_newton.upper_eta = 0.1
gauss_newton.save_intermediate = True
solution_low = mloc.solve()
logger.info("Solution: {}".format(solution_low.solution))
Starting solver MultiLevelIterativeSolver
Current option values:
abs_tol: 0.1
rel_tol: 0.1
max_iter: 10
Starting solver GaussNewton
Current option values:
abs_tol: 0.1
rel_tol: 0.1
max_iter: 20
Starting iteration: 0
Updating lower level variables...
Current residual ||J^T r||_2: 0.6658138241615206
Current ||r||_2: 0.6470057103540946
Current ||J||_2: 1.0317489703858296
Current cond(J): 1.0
Current allowed lower level tolerance: 0.02331910759868696
New x value:
[1.62546746]
New jac value:
[[ 0. ]
[-0.69271174]
[-0.64637097]
[-0.40848617]]
Starting iteration: 1
Updating lower level variables...
Current residual ||J^T r||_2: 0.09548903338269828
Current ||r||_2: 0.17109116894885706
Current ||J||_2: 0.5586616718068592
Current cond(J): 1.0
Current allowed lower level tolerance: 0.008452008146137857
Solution: [1.93142118]
[16]:
import tabulate
saved = solution_low.params['intermediates']
print(tabulate.tabulate(saved,headers='keys',tablefmt='latex_booktabs',floatfmt='.2e'))
\begin{tabular}{rrrrrrr}
\toprule
iteration & x & r & J & JTr & atol & cond \\
\midrule
0 & 1.00e+00 & 6.47e-01 & 1.03e+00 & 6.66e-01 & 1.00e-01 & 1.00e+00 \\
1 & 1.63e+00 & 1.71e-01 & 5.59e-01 & 9.55e-02 & 2.33e-02 & 1.00e+00 \\
& 1.93e+00 & & & & 8.45e-03 & \\
\bottomrule
\end{tabular}