<jupyter_start><jupyter_text>Mathematical Optimization for EngineersLab 11 - Full discretization methodIn this exercise, we want to solve the van der Pol oscillator optimal control problem using the full discretization approach. The problem with state constraint is as follows: $\begin{align}\min_{u(t)} \int_{t_{0}}^{t_{f}} x_1^2(t) &+ x_2^{2}(t) + u(t)^2 dt\\s.t.  \dot{x} &= f(x(t),u(t))\\x(t_{0}) &= x_{0} \\x_1(t) &\geq - 0.4 \quad \forall t \in [t_{0},t_{f}]\\-0.3 &\leq u(t) \leq 1.0\\\end{align}$For the full discretization appraoch, the differential equations are discretized usually by implicit Runge-Kutta methods. In this exercise, we will use the implicit Euler method due to its simplicity.  Imports<jupyter_code>import numpy as np

from scipy.optimize import minimize
from scipy.optimize import Bounds
from scipy.optimize import NonlinearConstraint

from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>Right-hand side function of ODE<jupyter_code>def eval_rhs(x, u):
    # your code here
    
    return f<jupyter_output><empty_output><jupyter_text>Constraints for optimization - discretized ODE<jupyter_code>def cons(y, x0, finalTime, M, n_u):
    # time step h (equidistant here)
    # your code here
    # number of states n_x
    # your code here
    ceq = np.zeros(n_x*M)

    # backward Euler y_(k+1) = y_(k) + h.f(t_(k+1),y_(k+1))
    
    # evalue RHS first time step f1
    # your code here
    
    # first element of constraint ceq[0:n_x]
    # your code here
        
    # int(z) for z in [] just converts all elements in the list to integer
    # it's an implementation detail because python only allows integer indices
    
    # second time step onward
    for i in # your code here. #Hint: np.arange
    
        # populate equality constraints using the backward Euler method
        # You could define indicies needed for the states and constraints
        # your code here
            
    return ceq<jupyter_output><empty_output><jupyter_text>objective<jupyter_code>def objective(y, x0):     
    # define objective
    # your code here
    return obj<jupyter_output><empty_output><jupyter_text>Setting up optimization problem and solve it<jupyter_code>n_x = 3 # dimension of state vector x 
n_u = 1 # dimension of control vector u

M=50 # Number of discretization intervals

x1_0 = 0 # initial value of x1
x2_0 = 1 # initial value of x2
x3_0 = 0 # initial value of x3

x_initial = [0.0, 1.0, 0.0]
finalTime = 5 # final time

x1_min = # your code here # lower bound on x1
u_min = # your code here # lower bound on u
u_max = # your code here # upper bound on u


n = # your code here # number of optimization variables
iu = (np.arange(3,n,4)).astype(int) # indices of control variables
ix1 = [int(z)-3 for z in iu] # indices of state x1
ix2 = [int(z)-2 for z in iu] # indices of state x2
ix3 = [int(z)-1 for z in iu] # indices of state x3

lb = -float("inf")*np.ones(n) # initialize lower bounds on all variables
ub = float("inf")*np.ones(n) # initialize upper bounds on all variables

lb[iu] = # your code here) # lower bound of control 
ub[iu] = # your code here # upper bound of control
lb[ix1] = # your code here # lower bound on x1



x_guess = np.zeros(n)  # initial guess of optimzation variable

x_guess[ix1] = x1_0  # initial guess for state x1
x_guess[ix2] = x2_0  # initial guess for state x2
x_guess[ix3] = x3_0  # initial guess for state x3

x_guess[iu] = 0.0  # initial guess for control u

bounds = Bounds(lb, ub)

consWithInitial = lambda y: cons(y, x_initial, finalTime, M, n_u)


nonlinear_constraints = NonlinearConstraint(consWithInitial, 0, 0)
# call optimization with bounds, constraints and objective
# your code here

print('Optimal function value is')
print(res.fun)<jupyter_output><empty_output><jupyter_text>Plot trajectories<jupyter_code>xopt = res.x

t = np.linspace(0, finalTime, M+1)
x1 = np.append([x1_0],list(xopt[np.arange(0,len(xopt),4)]))
x2 = np.append([x2_0],list(xopt[np.arange(1,len(xopt),4)]))
x3 = np.append([x3_0],list(xopt[np.arange(2,len(xopt),4)]))
u = np.append([0.0],list(xopt[np.arange(3,len(xopt),4)]))

plt.figure()
ax1 = plt.axes(xlim=(-4,4), ylim=(-4, 4))
plt.subplot(311)

plt.subplot(221)
plt.plot(t, x1)
plt.title('x1')

plt.subplot(222)
plt.plot(t, x2)
plt.title('x2')

plt.subplot(223)
plt.plot(t, x3)
plt.title('x3')
#
plt.subplot(224)
plt.plot(t, u)
plt.title('u')

plt.tight_layout()

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12092605079793212, 'fraction_numerical': 0.028320971004720162, 'mean_word_length': 4.103211009174312, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '81920', 'license': 'MIT', 'n_tokens_mistral': 1679, 'n_tokens_neox': 1617, 'n_words': 493, 'path': 'EngineeringOptimization/GitLab/Lab11.ipynb', 'repo_name': 'dalexa10/EngineeringDesignOptimization', 'size': '4449'}