<div>

</div>

<h1><center>Mathematical Optimization for Engineers</center></h1>
<h2><center>Lab 11 - Full discretization method</center></h2>

In 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


```python
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
```

### Right-hand side function of ODE


```python
def eval_rhs(x, u):
    # your code here
    
    return f
```

### Constraints for optimization - discretized ODE


```python
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
```

### objective


```python
def objective(y, x0):     
    # define objective
    # your code here
    return obj
```

### Setting up optimization problem and solve it


```python
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)
```

### Plot trajectories


```python
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()

```


```python

```

59939010e2d5259e437df3b8236beb000a41352a
ipynb
Jupyter Notebook
EngineeringOptimization/GitLab/Lab11.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
EngineeringOptimization/GitLab/Lab11.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
EngineeringOptimization/GitLab/Lab11.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
Qwen/Qwen-72B
1. YES
2. YES