<jupyter_start><jupyter_text>$m \ddot{x} + c \dot{x} + k x + sin(x) = u$ $\vec{x} = \begin{bmatrix}x \\\dot{x}\end{bmatrix}$ $\vec{u} = \begin{bmatrix} u\end{bmatrix}$ $\vec{y} = \vec{g}(\vec{x}) = \begin{bmatrix} x\end{bmatrix}$ $\ddot{x} = (-c \dot{x} - kx + u)/m$ $\dot{\vec{x}} = \vec{f}(\vec{x}) = \begin{bmatrix}\dot{x} \\(-c \dot{x} - kx - sin(x) + u)/m\end{bmatrix}$ $\dot{\vec{x}} = A \vec{x} + B \vec{u}$$\vec{y} = C \vec{x} + D \vec{u}$ $A = \dfrac{\partial \vec{f}}{\partial \vec{x}}$$B = \dfrac{\partial \vec{f}}{\partial \vec{u}}$$C = \dfrac{\partial \vec{g}}{\partial \vec{x}}$$D = \dfrac{\partial \vec{g}}{\partial \vec{u}}$<jupyter_code>m = ca.SX.sym('m')
c = ca.SX.sym('c')
k = ca.SX.sym('k')
p = ca.vertcat(m, c, k)

u = ca.SX.sym('u')
xv = ca.SX.sym('x', 2)
x = xv[0]
xd = xv[1]

y = x

xv_dot = ca.vertcat(xd, (-c*xd - k*x - ca.sin(x) + u + 3)/m)
xv_dot
f_rhs = ca.Function('rhs', [xv, u, p], [xv_dot], ['x', 'u', 'p'], ['x_dot'], {'jit': True})
f_rhs
f_rhs([1, 2], [0], [1, 2, 3])
import scipy.integrate
import numpy as np
tf = 10
res = scipy.integrate.solve_ivp(
    fun=lambda t, x: np.array(f_rhs(x, 0.0, [1, 2, 3])).reshape(-1),
    t_span=[0, tf],
    y0=[0, 0], t_eval=np.arange(0, tf, 0.1))
plt.plot(res['t'], res['y'][0, :]);
A = ca.jacobian(xv_dot, xv)
A
B = ca.jacobian(xv_dot, u)
B
C = ca.jacobian(y, xv)
C
D = ca.jacobian(y, u)
D
f_ss = ca.Function('f_ss', [xv, p], [A, B, C, D], ['x', 'p'], ['A', 'B', 'C', 'D'])
f_ss
import control
sys = control.ss(*f_ss([0, 0], [1, 2, 3]))
sys
f_rhs.generate('rhs.c')
#!cat rhs.c
s = control.TransferFunction([1, 0], [0, 1])
H = (s + 2)
control.rlocus(H*sys);
H*sys<jupyter_output><empty_output><jupyter_text>Linear Time Invariant Systems (LTI) * Transfer Functions: $G(s) = s/(s+1)$* State-space: $\dot{x} = Ax + Bu$, $y = Cx + Du$* Impulse response function: $g(t)$ * $\dot{x} = a_1 x + a_2 x + b u$, $y = c x + du$ Linear?  (Yes) Because A = A1 + A2* $\dot{x} = a_1 x + 3 + b u$, $y = c x + du$ Linear?  (No, not a linear system) * What u would balance this equation at x=0?  -> u0 = -3/b (trim input) For compensated dynamcis to be $G(s) = 1/(s+1)$,  u(x)=? * LTI $\implies$ zero in -> zero out $u(x) = (-a1 x - x - 3)/b$$\dot{x} = -x$  Trimming the MSD<jupyter_code>f_rhs([0, 0], [-3], [1, 2, 3])<jupyter_output><empty_output><jupyter_text>$\dot{x} = Ax + Bu$, $y = Cx + Du + 3$   (non-linear -> violates zero in zero out law) Trimming an aircraft means, finding where the rhs = 0,  or $f(t, x) = 0$, in order to do this we want to minimize$dot(f(t, x), f(t, x))$.<jupyter_code>def trim_function(xv_dot):
#    return xv_dot[0] + xv_dot[1] # BAD, will drive to -inf
    return xv_dot[0]**2 + xv_dot[1]**2<jupyter_output><empty_output><jupyter_text>This design problems find the state at which a given input will drive the sytem to.* x is the design vector* f is the objective function* p is a list of constant parameters* S is the solver itself<jupyter_code>nlp = {'x':xv, 'f':trim_function(xv_dot), 'p': ca.vertcat(p, u)}
S = ca.nlpsol('S', 'ipopt', nlp)
print(S)
S(x0=(0, 0), p=(1, 2, 3, 0))
nlp = {'x':u, 'f':trim_function(xv_dot), 'p': ca.vertcat(p, xv)}
S2 = ca.nlpsol('S', 'ipopt', nlp)
print(S2)
#res = S(x0=(0), p=(1, 2, 3, 0, 0))
#print('we need a trim input of {:f}'.format(float(res['x'])))

from casadi.tools.graph import dotgraph
from IPython.display import Image

def draw_graph(expr):
    return Image(dotgraph(expr).create_png())

S=2170.0; 
CBAR=17.5; 
MASS=5.0e3; 
IYY= 4.1e6;
TSTAT=6.0e4; 
DTDV = -38.0;
ZE = 2.0; 
CDCLS= 0.042;
CLA = 0.085; # per degree
CMA =-0.022; # per degree
CMDE =-.016; # per degree
CMQ =-16.0; # per radian
CMADOT= -6.0; # per radian
CLADOT= 0.0; # per radian
RTOD = 57.29578; 
GD=32.17;
#code 190
#4tran 725 standered atmospher model
import casadi as ca
THTL = ca.SX.sym('THTL')
ELEV = ca.SX.sym('ELEV')
XCG = ca.SX.sym('XCG')
LAND = ca.SX.sym('LAND')
PHI = ca.SX.sym('PHI')
u = ca.vertcat(THTL, ELEV, XCG, LAND, PHI)

VT = ca.SX.sym('VT')
ALPHA = ca.SX.sym('ALPHA')
THETA = ca.SX.sym('THETA')
Q = ca.SX.sym('Q')
H = ca.SX.sym('H')
PX = ca.SX.sym('PX')
RO = ca.SX(2.377e-3)
TFAC = ca.SX(1.0 - 0.703e-5 * H)
T = ca.SX(519*TFAC)
RHO = ca.SX(RO*(TFAC**4.14))#density
QBAR = ca.SX(0.5*RHO*VT*VT)#dynamic Pressure
CPS = ca.SX(1715.0*RHO*T)
#draw_graph(QBAR)
#no flaps
CLO = ca.SX(.20)
CDO = ca.SX(0.016)
CMO = ca.SX(0.05)
DCDG = ca.SX(0)
DCMG = ca.SX(0)
# Equations
QS = ca.SX(QBAR*S)
SALP = ca.sin(ALPHA+PHI)
CALP = ca.cos(ALPHA+PHI)
GAM = ca.SX(THETA-ALPHA)
SGAM = ca.sin(GAM)
CGAM = ca.cos(GAM)

THR = ca.SX((TSTAT+VT*DTDV) * ca.fmax(THTL,0))#Thrust
CL = ca.SX(CLO+CLA*ALPHA) #NONDIM. LIFT
CM = ca.SX(DCMG+CMO+CMA*ALPHA+CMDE*ELEV+CL* (XCG-.25)) #MOMENT
CD = ca.SX(DCDG+CDO+CDCLS*CL*CL) #DRAG POLAR
#EQU of motion
x1 = ((THR*CALP-QS*CD)/MASS - GD*SGAM)
x2 = ((-THR*SALP-QS*CL+MASS*(VT*Q+GD*CGAM))/(MASS*VT+QS*CLADOT))
x3 = (Q)
D = ((1/2)*CBAR*(CMQ*Q+CMADOT*x2/VT)) #PITCH DAMPING
x4 = ((QS*CBAR*(CM+D)+THR*ZE)/IYY) #Q-DOT
x5 = (VT*SGAM) #VERTICAL SPEED
x6 = (VT*CGAM) #ORIZNTal SPEED
totalbal=x1**2 + x2**2 + x4**2 + x5**2
trimm=ca.mmin(totalbal)
#draw_graph(trimm) 
x = ca.vertcat(x1, x2, x3, x4, x5, x6)
def objective(s, VT, H, GAMMA):
    THRTL = s[0]
    ELEV = s[1]
    ALPHA = s[2]
    PHI = s[3]
    
    POS = 0
    Q = 0
    XCG = 0.25
    THETA = ALPHA+GAMMA
    
    x = ca.vertcat(thtl, elev, xcg, land)
    
    x = ca.vertcat
nlp = {'x': s, 'f': objectives(s, VT=500, H=0, GAMMA=0)}
S =<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.23884707178263265, 'fraction_numerical': 0.04020561777125023, 'mean_word_length': 3.8599464763603923, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '135050', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2861, 'n_tokens_neox': 2712, 'n_words': 622, 'path': 'lectures/4-Casadi-MSD-Copy1.ipynb', 'repo_name': 'HenryTingle/aae497-f19', 'size': '5447'}