<div>

</div>
<h1><center>Mathematical Optimization for Engineers</center></h1>
<h2><center>Lab 12 - Direct single shooting</center></h2>

In this exercise, we want to solve the van der Pol oscillator optimal control problem using the direct single shooting approach. 
The problem without state constraints 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} \\
u_{lb} \leq u(t) \leq u_{ub}\\
\end{align}
$

For direct single shooting, the differential equations are solved by an integration scheme, which is needed to be set-up. 



```python
import numpy as np

from scipy.integrate import odeint, solve_ivp

from scipy.optimize import fmin,fmin_bfgs, Bounds, minimize, NonlinearConstraint

from matplotlib import pyplot as plt
from matplotlib import animation
```


```python
def eval_rhs(t, x, u, p):
    
    # your code

    
    return f
```


```python
def objective_ode_integration(p, func, x0, finalTime, output_state,
                              is_path_constraint, doPlot):
    # inputs are as follows
    # p = parameter vector of control discretization
    # func: right-hand side of ODE f(t, x, u, p)
    # x0 = initial condition
    # finalTime = tN
    # doPlot = boolean variable, if true, plot solution   

    N = len(p) # Number of discretization intervals
    h = finalTime /N # length of a single interval
    n_x = len(x0)  # number of states
    x0_cur = x0
    y_all = np.array([x0]).transpose() 
    t_all = np.array([0])
    output = list()

    for i in range(N): # loop over discretization intervals
        
        # solve initial value problem
        # your code
        
        # update initial conditions
        # your code

        y_all = np.concatenate((y_all, sol.y), axis=1)
        t_all = np.concatenate((t_all, sol.t), axis=0)
        output.append(y[output_state, -1])



    if doPlot:
        y_all = np.array(y_all)
        plt.subplot(221)
        plt.plot(t_all, y_all[0, :])
        plt.title('x1')

        plt.subplot(222)
        plt.plot(t_all, y_all[1, :])
        plt.title('x2')

        plt.subplot(223)
        plt.plot(t_all, y_all[2, :])
        plt.title('x3')
        #
        plt.subplot(224)
        plt.plot(np.linspace(0, finalTime, N), p)
        plt.title('u')

        plt.tight_layout()
        plt.show()
        
    if not is_path_constraint: 
        output = output[-1]
    else:
        # constraint hard coded here to allow for 
        # constraint definition with arguments
        # in scipy, c(x) >= 0
        # your code here

    return output
```


```python
objective_ode_integration([1.0, 0.5, 0.2, 0.4], eval_rhs, [0, 1, 0], 5, [], [], True)
```


```python
finalTime = 5
N = 25 # Number of discretization intervals for the control
doPlot = False
p0 = 1.0*np.ones(N) # initial guess for optimization
x0 = [0, 1, 0]

bnds = Bounds(-0.3, 1.0)

# define objective
output_state = 2  # for objective
is_path_constraint = False
# call optimization with arguments 
# 'eval_rhs, x0, finalTime, output_state, is_path_constraint, doPlot'
# your code here

# call ODE integration to plot trajectories with optimal solution of controls
# your code here

```

## Adding state path constraint
After solvindg the problem without state constraints, we want to include the path constraint on the state $x_1$:
$ \begin{align}
x_1(t) \geq - 0.4 \forall t \in [t_{0},t_{f}]
\end{align}$
For this task, we must add this constraint to the optimization problem. One way is to define the constraint function as constraint dictionary in Python. 


```python
# in scipy, the inequality constraints are defined as
# c(x) >= 0 
cons = ({'type': 'ineq',
       'fun': objective_ode_integration,
       'args': (eval_rhs, x0, finalTime, 0, True, False)       
       })

output_state = 2  # for objective
is_path_constraint = False

# call optimization with arguments 
# 'eval_rhs, x0, finalTime, output_state, is_path_constraint, doPlot'
# your code here


# call ODE integration to plot trajectories with optimal solution of controls
# your code here
```

## What could be made more efficient for this optimization? 


```python

```

77789a7f3703896954daa6f68703b2c462b821d0
ipynb
Jupyter Notebook
EngineeringOptimization/GitLab/Lab12.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
EngineeringOptimization/GitLab/Lab12.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
EngineeringOptimization/GitLab/Lab12.ipynb
dalexa10/EngineeringDesignOptimization
eb5b5e4edd773aef629f59aea8a9771af41bd224
MIT
Qwen/Qwen-72B
1. YES
2. YES