# Trimming

Trimming, or determining an operating point, of a dynamic system, involves finding a state and input to the system dynamics that satisfy user-defined constraints. Trimming is commonly used as a first step in determining a point about which to linearize a nonlinear system.

A tutorial on trimming in the form of a webinar is available below

## Trim vs Equilibrium Points:

A special case of trimming is finding an equilibrium point of a dynamic system. An equilibrium point is a feasible steady state of the dynamic system, i.e., a state at which all the derivatives are zero, $f(x, u, p, t) = 0$. Trimming can be applied to many dynamic system formalisms; this document will discuss ODE systems and DAE systems.

## ODE Systems

Consider an ODE system $\dot x = f(x,u,p,t)$ where $x \in \mathbb{R}^n$ are the state variables, $u \in \mathbb{R}^m$ control input variables, $p$ are the parameters and $f:\mathbb{R}^{n} \times \mathbb{R}^{m} \mapsto \mathbb{R}^n$ is a potentially nonlinear function. An equilibrium point $(x_{eq}, u_{eq})$ is a solution to the set of nonlinear equations $f(x_{eq}, u_{eq}, p, t) = \mathbf{0}$. This point is generally not unique since the equation system has $n+m$ variables and $n$ equations.

Various methods may be used to determine the equilibrium point: For instance, one may use the a root-finding algorithm to solve the system of nonlinear equations above. In addition, a constrained-optimization problem may be formulated which in general is a nonconvex nonlinear program (NLP). Alternatively, an unconstrained optimization problem may be formulated by using a penalty function to dualize the constraints.

Trim points generalize the concept of equilibrium points: A user may want to impose additional hard or soft constraints (defined by equations or inequalities involving the states of the system). In addition, a user may specify certain derivatives to be non-zero. Unlike equilibrium points which may be determined using root-finding algorithms, it is usually necessary to formulate a constrained or unconstrained optimization problem to determine trim points as is done in the developed trim function.

## Trimming Workflow:

1. The dynamic system here named sys is defined symbolically using the ModelingToolkit framework, thus sys is of type AbstractSystem. This system should not be simplified by the user by calling structural_simplify as this is done inside the trim function.
• Next the trim function is called. If the system is unbalanced, i.e., if the number of states exceeds the number of equations as a result of having certain control inputs, it is necessary to specify an array of control inputs and outputs using the inputs and outputs arguments.
• Optionally, to determine a trim point, the user may specify certain desired states, hard or soft inequality or equality constraints using the corresponding arguments. If none of these are specified, the function returns an equilibrium point.
• By default, the internals of the trim function involve formulating a constrained optimization problem. The user may want to specify tailored penalty multipliers or use the penalty function approach to dualize the contraints such that an unconstrained optimization problem is formulated instead. In either case, it is essential to specify an appropriate optimization solver.

## Example: The nonlinear Research Civil Aircraft Model (RCAM)

The trim functionality is illustrated using the nonlinear Research Civil Aircraft Model (RCAM) developed by GARTEUR and using the derivations and errata presented by Prof. Christopher Lum. An open-source implementation of this model is also presented in PSim-RCAM.

Figure 1: Aircraft axes and control surfaces. Figure attained from Wikimedia Commons with CC by 4.0 license,

The RCAM model has the following 9 states:

• $$$u, v, w$$$: the translational velocity [m/s] in the body frame in the x, y and z-axes respectively,
• $$$p, q, r$$$: the rotational velocity [rad/s] in the body frame about the x (roll/bank), y (pitch) and z (yaw)-axes respectively,
• $$$ϕ, θ, ψ$$$: the rotation angles [rad] in the body frame about the x (roll/bank), y (pitch) and z (yaw)-axes respectively,

and the following 5 control inputs:

• $$$u_A$$$: the aileron angle [rad]
• $$$u_T$$$: the tail angle [rad]
• $$$u_R$$$: the rudder angle [rad]
• $$$u_{E1}$$$: the left engine/throttle position [rad]
• $$$u_{E2}$$$: the right engine/throttle position [rad]

The mechanistic model is defined symbolically as an ODESystem in ModelingToolkit using the 9-step procedure and parameters detailed in this lecture.

1.a. First, the parameters for the model are defined as follows:

using ModelingToolkit, OrdinaryDiffEq, ModelingToolkit.IfElse, JuliaSimControl
using LinearAlgebra:cross, inv
using Plots

# Define a dictionary of constants

rcam_constants = Dict(:m => 120000, :c̄ => 6.6, :lt => 24.8, :S => 260, :St => 64, :Xcg => 0.23 * 6.6,:Ycg => 0.0, :Zcg => 0.10 * 6.6, :Xac => 0.12 * 6.6, :Yac => 0.0, :Zac => 0.0, :Xapt1 => 0, :Yapt1 => -7.94, :Zapt1 => -1.9, :Xapt2 => 0, :Yapt2 => 7.94, :Zapt2 => -1.9, :g => 9.81, :depsda => 0.25, :α_L0 => deg2rad(-11.5), :n => 5.5, :a3 => -768.5, :a2 => 609.2, :a1 => -155.2, :a0 => 15.212, :α_switch => deg2rad(14.5), :ρ => 1.225)

Ib_mat =  [40.07          0.0         -2.0923
0.0            64.0        0.0
-2.0923        0.0         99.92]*rcam_constants[:m]

Ib = Ib_mat
invIb = inv(Ib_mat)

ps = @parameters(
ρ             = rcam_constants[:ρ], [description="kg/m3 - air density"],
m             = rcam_constants[:m], [description="kg - total mass"],
c̄          = rcam_constants[:c̄], [description="m - mean aerodynamic chord"],
lt            = rcam_constants[:lt], [description="m - tail AC distance to CG"],
S             = rcam_constants[:S], [description="m2 - wing area"],
St            = rcam_constants[:St], [description="m2 - tail area"],
Xcg           = rcam_constants[:Xcg], [description="m - x pos of CG in Fm"],
Ycg           = rcam_constants[:Ycg], [description="m - y pos of CG in Fm"],
Zcg           = rcam_constants[:Zcg], [description="m - z pos of CG in Fm"],
Xac           = rcam_constants[:Xac], [description="m - x pos of aerodynamic center in Fm"],
Yac           = rcam_constants[:Yac], [description="m - y pos of aerodynamic center in Fm"],
Zac           = rcam_constants[:Zac], [description="m - z pos of aerodynamic center in Fm"],
Xapt1         = rcam_constants[:Xapt1], [description="m - x position of engine 1 in Fm"],
Yapt1         = rcam_constants[:Yapt1], [description="m - y position of engine 1 in Fm"],
Zapt1         = rcam_constants[:Zapt1], [description="m - z position of engine 1 in Fm"],
Xapt2         = rcam_constants[:Xapt2], [description="m - x position of engine 2 in Fm"],
Yapt2         = rcam_constants[:Yapt2], [description="m - y position of engine 2 in Fm"],
Zapt2         = rcam_constants[:Zapt2], [description="m - z position of engine 2 in Fm"],
g             = rcam_constants[:g], [description="m/s2 - gravity"],
α_L0      = rcam_constants[:α_L0], [description="rad - zero lift AOA"],
n             = rcam_constants[:n], [description="adm - slope of linear ragion of lift slope"],
a3            = rcam_constants[:a3], [description="adm - coeff of α^3"],
a2            = rcam_constants[:a2], [description="adm -  - coeff of α^2"],
a1            = rcam_constants[:a1], [description="adm -  - coeff of α^1"],
a0            = rcam_constants[:a0], [description="adm -  - coeff of α^0"],
α_switch  = rcam_constants[:α_switch], [description="rad - kink point of lift slope"]
)

1.b. Next the 9 states and 5 control variables are defined:

@parameters t
D =  Differential(t)

# States
V_b = @variables(
u(t), [description="translational velocity along x-axis [m/s]"],
v(t), [description="translational velocity along y-axis [m/s]"],
w(t), [description="translational velocity along z-axis [m/s]"]
)

wbe_b = @variables(
)

rot = @variables(
)

# Controls
U = @variables(
)

1.c. Next, the auxiliary variables and intermediate system of equations is defined:

# Auxiliary Variables to define model.

Auxiliary_vars = @variables Va(t) α(t) β(t) Q(t) CL_wb(t) ϵ(t) α_t(t) CL_t(t) CL(t) CD(t) CY(t) F1(t) F2(t)

@variables FA_s(t)[1:3] C_bs(t)[1:3,1:3] FA_b(t)[1:3] eta(t)[1:3] dCMdx(t)[1:3, 1:3] dCMdu(t)[1:3, 1:3] CMac_b(t)[1:3] MAac_b(t)[1:3] rcg_b(t)[1:3] rac_b(t)[1:3] MAcg_b(t)[1:3] FE1_b(t)[1:3] FE2_b(t)[1:3] FE_b(t)[1:3] mew1(t)[1:3] mew2(t)[1:3] MEcg1_b(t)[1:3] MEcg2_b(t)[1:3] MEcg_b(t)[1:3] g_b(t)[1:3] Fg_b(t)[1:3] F_b(t)[1:3] Mcg_b(t)[1:3] H_phi(t)[1:3,1:3]

# Scalarizing all the array variables.

FA_s = collect(FA_s)
C_bs = collect(C_bs)
FA_b = collect(FA_b)
eta = collect(eta)
dCMdx = collect(dCMdx)
dCMdu = collect(dCMdu)
CMac_b =  collect(CMac_b)
MAac_b = collect(MAac_b)
rcg_b = collect(rcg_b)
rac_b = collect(rac_b)
MAcg_b = collect(MAcg_b)
FE1_b = collect(FE1_b)
FE2_b = collect(FE2_b)
FE_b = collect(FE_b)
mew1 = collect(mew1)
mew2 = collect(mew2)
MEcg1_b = collect(MEcg1_b)
MEcg2_b = collect(MEcg2_b)
MEcg_b = collect(MEcg_b)
g_b = collect(g_b)
Fg_b = collect(Fg_b)
F_b = collect(F_b)
Mcg_b = collect(Mcg_b)
H_phi = collect(H_phi)

array_vars = vcat(vec(FA_s), vec(C_bs), vec(FA_b), vec(eta), vec(dCMdx), vec(dCMdu), vec(CMac_b), vec(MAac_b), vec(rcg_b), vec(rac_b), vec(MAcg_b), vec(FE1_b), vec(FE2_b), vec(FE_b), vec(mew1), vec(mew2), vec(MEcg1_b), vec(MEcg2_b), vec(MEcg_b), vec(g_b), vec(Fg_b), vec(F_b), vec(Mcg_b), vec(H_phi))

eqns =[
# Step 1. Intermediate variables
# Airspeed
Va ~ sqrt(u^2 + v^2 + w^2)

# α and β
α ~ atan(w,u)
β ~ asin(v/Va)

# dynamic pressure
Q ~ 0.5*ρ*Va^2

# Step 2. Aerodynamic Force Coefficients
# CL - wing + body
#CL_wb ~  n*(α - α_L0)

CL_wb ~ IfElse.ifelse(α <= α_switch, n*(α - α_L0), a3*α^3 + a2*α^2 + a1*α + a0)

# CL thrust
ϵ ~ depsda*(α - α_L0)
α_t ~ α - ϵ + uT + 1.3*q*lt/Va
CL_t ~ 3.1*(St/S) * α_t

# Total CL
CL ~ CL_wb + CL_t

# Total CD
CD ~ 0.13 + 0.07 * (n*α + 0.654)^2

# Total CY
CY ~ -1.6*β + 0.24*uR

# Step 3. Dimensional Aerodynamic Forces
# Forces in F_s
FA_s .~ [-CD * Q * S
CY * Q * S
-CL * Q * S]

# rotate forces to body axis (F_b)
vec(C_bs .~ [cos(α)      0.0      -sin(α)
0.0             1.0      0.0
sin(α)      0.0      cos(α)])

FA_b .~ C_bs*FA_s

# Step 4. Aerodynamic moment coefficients about AC
# moments in F_b
eta .~ [ -1.4 * β
-0.59 - (3.1 * (St * lt) / (S * c̄)) * (α - ϵ)
(1 - α * (180 / (15 * π))) * β
]

vec(dCMdx .~ (c̄ / Va)*        [-11.0              0.0                           5.0
0.0   (-4.03 * (St * lt^2) / (S * c̄^2))        0.0
1.7                 0.0                          -11.5])

vec(dCMdu .~ [-0.6                   0.0                 0.22
0.0   (-3.1 * (St * lt) / (S * c̄))      0.0
0.0                    0.0                -0.63])

# CM about AC in Fb
CMac_b .~ eta + dCMdx*wbe_b + dCMdu*[uA
uT
uR]

# Step 5. Aerodynamic moment about AC
# normalize to aerodynamic moment
MAac_b .~ CMac_b * Q * S * c̄

# Step 6. Aerodynamic moment about CG
rcg_b .~    [Xcg
Ycg
Zcg]

rac_b .~ [Xac
Yac
Zac]

MAcg_b .~ MAac_b + cross(FA_b, rcg_b - rac_b)

# Step 7. Engine force and moment
# thrust
F1 ~ uE_1 * m * g
F2 ~ uE_2 * m * g

# thrust vectors (assuming aligned with x axis)
FE1_b .~ [F1
0
0]

FE2_b .~ [F2
0
0]

FE_b .~ FE1_b + FE2_b

# engine moments
mew1 .~  [Xcg - Xapt1
Yapt1 - Ycg
Zcg - Zapt1]

mew2 .~ [ Xcg - Xapt2
Yapt2 - Ycg
Zcg - Zapt2]

MEcg1_b .~ cross(mew1, FE1_b)
MEcg2_b .~ cross(mew2, FE2_b)

MEcg_b .~ MEcg1_b + MEcg2_b

# Step 8. Gravity effects
g_b .~ [-g * sin(θ)
g * cos(θ) * sin(ϕ)
g * cos(θ) * cos(ϕ)]

Fg_b .~ m * g_b

# Step 9: State derivatives

# form F_b and calculate u, v, w dot
F_b .~ Fg_b + FE_b + FA_b

D.(V_b) .~ (1 / m)*F_b - cross(wbe_b, V_b)

# form Mcg_b and calc p, q r dot
Mcg_b .~ MAcg_b + MEcg_b

D.(wbe_b) .~ invIb*(Mcg_b - cross(wbe_b, Ib*wbe_b))

#phi, theta, psi dot
vec(H_phi .~ [1.0         sin(ϕ)*tan(θ)       cos(ϕ)*tan(θ)
0.0         cos(ϕ)              -sin(ϕ)
0.0         sin(ϕ)/cos(θ)       cos(ϕ)/cos(θ)])

D.(rot) .~  H_phi*wbe_b
]

1.d. Lastly, an ODESystem <: AbstractSystem is formed. The gives a symbolic ODE system with 118 equations and 123 variables. Thus, there are 5 control inputs to the ODE system.

all_vars = vcat(V_b,wbe_b,rot,U, Auxiliary_vars, array_vars)
@named rcam_model = ODESystem(eqns, t, all_vars, ps)
\begin{align} \mathrm{Va}\left( t \right) =& \sqrt{\left( u\left( t \right) \right)^{2} + \left( v\left( t \right) \right)^{2} + \left( w\left( t \right) \right)^{2}} \\ \alpha\left( t \right) =& \arctan\left( w\left( t \right), u\left( t \right) \right) \\ \beta\left( t \right) =& \arcsin\left( \frac{v\left( t \right)}{\mathrm{Va}\left( t \right)} \right) \\ Q\left( t \right) =& 0.5 \left( \mathrm{Va}\left( t \right) \right)^{2} \rho \\ \mathrm{CL}_{wb}\left( t \right) =& \mathrm{ifelse}\left( \left( \alpha\left( t \right) \leq \alpha_{switch} \right), n \left( - \alpha_{L0} + \alpha\left( t \right) \right), a0 + \left( \alpha\left( t \right) \right)^{2} a2 + \left( \alpha\left( t \right) \right)^{3} a3 + a1 \alpha\left( t \right) \right) \\ \epsilon\left( t \right) =& depsda \left( - \alpha_{L0} + \alpha\left( t \right) \right) \\ \alpha_{t}\left( t \right) =& \frac{1.3 lt q\left( t \right)}{\mathrm{Va}\left( t \right)} - \epsilon\left( t \right) + \mathrm{uT}\left( t \right) + \alpha\left( t \right) \\ \mathrm{CL}_{t}\left( t \right) =& \frac{3.1 St \alpha_{t}\left( t \right)}{S} \\ \mathrm{CL}\left( t \right) =& \mathrm{CL}_{t}\left( t \right) + \mathrm{CL}_{wb}\left( t \right) \\ \mathrm{CD}\left( t \right) =& 0.13 + 0.07 \left( 0.654 + n \alpha\left( t \right) \right)^{2} \\ \mathrm{CY}\left( t \right) =& 0.24 \mathrm{uR}\left( t \right) - 1.6 \beta\left( t \right) \\ FA_{s(t)_1} =& - S \mathrm{CD}\left( t \right) Q\left( t \right) \\ FA_{s(t)_2} =& S \mathrm{CY}\left( t \right) Q\left( t \right) \\ FA_{s(t)_3} =& - S \mathrm{CL}\left( t \right) Q\left( t \right) \\ C_{bs(t)_{1}ˏ_1} =& \cos\left( \alpha\left( t \right) \right) \\ C_{bs(t)_{2}ˏ_1} =& 0.0 \\ C_{bs(t)_{3}ˏ_1} =& \sin\left( \alpha\left( t \right) \right) \\ C_{bs(t)_{1}ˏ_2} =& 0.0 \\ C_{bs(t)_{2}ˏ_2} =& 1.0 \\ C_{bs(t)_{3}ˏ_2} =& 0.0 \\ C_{bs(t)_{1}ˏ_3} =& - \sin\left( \alpha\left( t \right) \right) \\ C_{bs(t)_{2}ˏ_3} =& 0.0 \\ C_{bs(t)_{3}ˏ_3} =& \cos\left( \alpha\left( t \right) \right) \\ FA_{b(t)_1} =& FA_{s(t)_1} C_{bs(t)_{1}ˏ_1} + FA_{s(t)_2} C_{bs(t)_{1}ˏ_2} + FA_{s(t)_3} C_{bs(t)_{1}ˏ_3} \\ FA_{b(t)_2} =& FA_{s(t)_1} C_{bs(t)_{2}ˏ_1} + FA_{s(t)_2} C_{bs(t)_{2}ˏ_2} + FA_{s(t)_3} C_{bs(t)_{2}ˏ_3} \\ FA_{b(t)_3} =& FA_{s(t)_1} C_{bs(t)_{3}ˏ_1} + FA_{s(t)_2} C_{bs(t)_{3}ˏ_2} + FA_{s(t)_3} C_{bs(t)_{3}ˏ_3} \\ eta(t)_1 =& - 1.4 \beta\left( t \right) \\ eta(t)_2 =& -0.59 + \frac{ - 3.1 St lt \left( - \epsilon\left( t \right) + \alpha\left( t \right) \right)}{S \textnormal{\={c}}} \\ eta(t)_3 =& \left( 1 - 3.8197186342054885 \alpha\left( t \right) \right) \beta\left( t \right) \\ dCMdx(t)_{1}ˏ_1 =& \frac{ - 11.0 \textnormal{\={c}}}{\mathrm{Va}\left( t \right)} \\ dCMdx(t)_{2}ˏ_1 =& 0.0 \\ dCMdx(t)_{3}ˏ_1 =& \frac{1.7 \textnormal{\={c}}}{\mathrm{Va}\left( t \right)} \\ dCMdx(t)_{1}ˏ_2 =& 0.0 \\ dCMdx(t)_{2}ˏ_2 =& \frac{ - 4.03 lt^{2} St}{S \textnormal{\={c}} \mathrm{Va}\left( t \right)} \\ dCMdx(t)_{3}ˏ_2 =& 0.0 \\ dCMdx(t)_{1}ˏ_3 =& \frac{5.0 \textnormal{\={c}}}{\mathrm{Va}\left( t \right)} \\ dCMdx(t)_{2}ˏ_3 =& 0.0 \\ dCMdx(t)_{3}ˏ_3 =& \frac{ - 11.5 \textnormal{\={c}}}{\mathrm{Va}\left( t \right)} \\ dCMdu(t)_{1}ˏ_1 =& -0.6 \\ dCMdu(t)_{2}ˏ_1 =& 0.0 \\ dCMdu(t)_{3}ˏ_1 =& 0.0 \\ dCMdu(t)_{1}ˏ_2 =& 0.0 \\ dCMdu(t)_{2}ˏ_2 =& \frac{ - 3.1 St lt}{S \textnormal{\={c}}} \\ dCMdu(t)_{3}ˏ_2 =& 0.0 \\ dCMdu(t)_{1}ˏ_3 =& 0.22 \\ dCMdu(t)_{2}ˏ_3 =& 0.0 \\ dCMdu(t)_{3}ˏ_3 =& -0.63 \\ CMac_{b(t)_1} =& dCMdu(t)_{1}ˏ_1 \mathrm{uA}\left( t \right) + dCMdu(t)_{1}ˏ_2 \mathrm{uT}\left( t \right) + dCMdx(t)_{1}ˏ_1 p\left( t \right) + dCMdx(t)_{1}ˏ_2 q\left( t \right) + dCMdx(t)_{1}ˏ_3 r\left( t \right) + dCMdu(t)_{1}ˏ_3 \mathrm{uR}\left( t \right) + eta(t)_1 \\ CMac_{b(t)_2} =& dCMdx(t)_{2}ˏ_1 p\left( t \right) + dCMdu(t)_{2}ˏ_1 \mathrm{uA}\left( t \right) + dCMdx(t)_{2}ˏ_2 q\left( t \right) + dCMdu(t)_{2}ˏ_2 \mathrm{uT}\left( t \right) + dCMdu(t)_{2}ˏ_3 \mathrm{uR}\left( t \right) + dCMdx(t)_{2}ˏ_3 r\left( t \right) + eta(t)_2 \\ CMac_{b(t)_3} =& dCMdx(t)_{3}ˏ_1 p\left( t \right) + dCMdx(t)_{3}ˏ_2 q\left( t \right) + dCMdx(t)_{3}ˏ_3 r\left( t \right) + dCMdu(t)_{3}ˏ_1 \mathrm{uA}\left( t \right) + dCMdu(t)_{3}ˏ_2 \mathrm{uT}\left( t \right) + dCMdu(t)_{3}ˏ_3 \mathrm{uR}\left( t \right) + eta(t)_3 \\ MAac_{b(t)_1} =& S \textnormal{\={c}} Q\left( t \right) CMac_{b(t)_1} \\ MAac_{b(t)_2} =& S \textnormal{\={c}} Q\left( t \right) CMac_{b(t)_2} \\ MAac_{b(t)_3} =& S \textnormal{\={c}} Q\left( t \right) CMac_{b(t)_3} \\ rcg_{b(t)_1} =& Xcg \\ rcg_{b(t)_2} =& Ycg \\ rcg_{b(t)_3} =& Zcg \\ rac_{b(t)_1} =& Xac \\ rac_{b(t)_2} =& Yac \\ rac_{b(t)_3} =& Zac \\ MAcg_{b(t)_1} =& \left( - rac_{b(t)_3} + rcg_{b(t)_3} \right) FA_{b(t)_2} + \left( - rcg_{b(t)_2} + rac_{b(t)_2} \right) FA_{b(t)_3} + MAac_{b(t)_1} \\ MAcg_{b(t)_2} =& \left( - rcg_{b(t)_3} + rac_{b(t)_3} \right) FA_{b(t)_1} + \left( - rac_{b(t)_1} + rcg_{b(t)_1} \right) FA_{b(t)_3} + MAac_{b(t)_2} \\ MAcg_{b(t)_3} =& \left( - rcg_{b(t)_1} + rac_{b(t)_1} \right) FA_{b(t)_2} + \left( - rac_{b(t)_2} + rcg_{b(t)_2} \right) FA_{b(t)_1} + MAac_{b(t)_3} \\ \mathrm{F1}\left( t \right) =& g m \mathrm{uE}_{1}\left( t \right) \\ \mathrm{F2}\left( t \right) =& g m \mathrm{uE}_{2}\left( t \right) \\ FE1_{b(t)_1} =& \mathrm{F1}\left( t \right) \\ FE1_{b(t)_2} =& 0 \\ FE1_{b(t)_3} =& 0 \\ FE2_{b(t)_1} =& \mathrm{F2}\left( t \right) \\ FE2_{b(t)_2} =& 0 \\ FE2_{b(t)_3} =& 0 \\ FE_{b(t)_1} =& FE1_{b(t)_1} + FE2_{b(t)_1} \\ FE_{b(t)_2} =& FE1_{b(t)_2} + FE2_{b(t)_2} \\ FE_{b(t)_3} =& FE1_{b(t)_3} + FE2_{b(t)_3} \\ mew1(t)_1 =& Xcg - Xapt1 \\ mew1(t)_2 =& Yapt1 - Ycg \\ mew1(t)_3 =& Zcg - Zapt1 \\ mew2(t)_1 =& Xcg - Xapt2 \\ mew2(t)_2 =& Yapt2 - Ycg \\ mew2(t)_3 =& Zcg - Zapt2 \\ MEcg1_{b(t)_1} =& FE1_{b(t)_3} mew1(t)_2 - FE1_{b(t)_2} mew1(t)_3 \\ MEcg1_{b(t)_2} =& FE1_{b(t)_1} mew1(t)_3 - FE1_{b(t)_3} mew1(t)_1 \\ MEcg1_{b(t)_3} =& FE1_{b(t)_2} mew1(t)_1 - FE1_{b(t)_1} mew1(t)_2 \\ MEcg2_{b(t)_1} =& FE2_{b(t)_3} mew2(t)_2 - FE2_{b(t)_2} mew2(t)_3 \\ MEcg2_{b(t)_2} =& FE2_{b(t)_1} mew2(t)_3 - FE2_{b(t)_3} mew2(t)_1 \\ MEcg2_{b(t)_3} =& FE2_{b(t)_2} mew2(t)_1 - FE2_{b(t)_1} mew2(t)_2 \\ MEcg_{b(t)_1} =& MEcg1_{b(t)_1} + MEcg2_{b(t)_1} \\ MEcg_{b(t)_2} =& MEcg1_{b(t)_2} + MEcg2_{b(t)_2} \\ MEcg_{b(t)_3} =& MEcg1_{b(t)_3} + MEcg2_{b(t)_3} \\ g_{b(t)_1} =& - g \sin\left( \theta\left( t \right) \right) \\ g_{b(t)_2} =& g \cos\left( \theta\left( t \right) \right) \sin\left( \phi\left( t \right) \right) \\ g_{b(t)_3} =& g \cos\left( \theta\left( t \right) \right) \cos\left( \phi\left( t \right) \right) \\ Fg_{b(t)_1} =& m g_{b(t)_1} \\ Fg_{b(t)_2} =& m g_{b(t)_2} \\ Fg_{b(t)_3} =& m g_{b(t)_3} \\ F_{b(t)_1} =& FA_{b(t)_1} + FE_{b(t)_1} + Fg_{b(t)_1} \\ F_{b(t)_2} =& FA_{b(t)_2} + FE_{b(t)_2} + Fg_{b(t)_2} \\ F_{b(t)_3} =& FA_{b(t)_3} + FE_{b(t)_3} + Fg_{b(t)_3} \\ \frac{\mathrm{d} u\left( t \right)}{\mathrm{d}t} =& r\left( t \right) v\left( t \right) + \frac{F_{b(t)_1}}{m} - q\left( t \right) w\left( t \right) \\ \frac{\mathrm{d} v\left( t \right)}{\mathrm{d}t} =& p\left( t \right) w\left( t \right) + \frac{F_{b(t)_2}}{m} - r\left( t \right) u\left( t \right) \\ \frac{\mathrm{d} w\left( t \right)}{\mathrm{d}t} =& q\left( t \right) u\left( t \right) + \frac{F_{b(t)_3}}{m} - p\left( t \right) v\left( t \right) \\ Mcg_{b(t)_1} =& MAcg_{b(t)_1} + MEcg_{b(t)_1} \\ Mcg_{b(t)_2} =& MAcg_{b(t)_2} + MEcg_{b(t)_2} \\ Mcg_{b(t)_3} =& MAcg_{b(t)_3} + MEcg_{b(t)_3} \\ \frac{\mathrm{d} p\left( t \right)}{\mathrm{d}t} =& 4.3595940921659655e-9 \left( \left( 4.8084e6 p\left( t \right) - 251075.99999999997 r\left( t \right) \right) q\left( t \right) - 7.68e6 p\left( t \right) q\left( t \right) + Mcg_{b(t)_3} \right) + 2.0819702800230527e-7 \left( \left( 251075.99999999997 p\left( t \right) - 1.19904e7 r\left( t \right) \right) q\left( t \right) + 7.68e6 q\left( t \right) r\left( t \right) + Mcg_{b(t)_1} \right) \\ \frac{\mathrm{d} q\left( t \right)}{\mathrm{d}t} =& 1.3020833333333334e-7 \left( \left( - 251075.99999999997 p\left( t \right) + 1.19904e7 r\left( t \right) \right) p\left( t \right) + \left( - 4.8084e6 p\left( t \right) + 251075.99999999997 r\left( t \right) \right) r\left( t \right) + Mcg_{b(t)_2} \right) \\ \frac{\mathrm{d} r\left( t \right)}{\mathrm{d}t} =& 8.349134219427914e-8 \left( \left( 4.8084e6 p\left( t \right) - 251075.99999999997 r\left( t \right) \right) q\left( t \right) - 7.68e6 p\left( t \right) q\left( t \right) + Mcg_{b(t)_3} \right) + 4.3595940921659655e-9 \left( \left( 251075.99999999997 p\left( t \right) - 1.19904e7 r\left( t \right) \right) q\left( t \right) + 7.68e6 q\left( t \right) r\left( t \right) + Mcg_{b(t)_1} \right) \\ H_{phi(t)_{1}ˏ_1} =& 1.0 \\ H_{phi(t)_{2}ˏ_1} =& 0.0 \\ H_{phi(t)_{3}ˏ_1} =& 0.0 \\ H_{phi(t)_{1}ˏ_2} =& \tan\left( \theta\left( t \right) \right) \sin\left( \phi\left( t \right) \right) \\ H_{phi(t)_{2}ˏ_2} =& \cos\left( \phi\left( t \right) \right) \\ H_{phi(t)_{3}ˏ_2} =& \frac{\sin\left( \phi\left( t \right) \right)}{\cos\left( \theta\left( t \right) \right)} \\ H_{phi(t)_{1}ˏ_3} =& \tan\left( \theta\left( t \right) \right) \cos\left( \phi\left( t \right) \right) \\ H_{phi(t)_{2}ˏ_3} =& - \sin\left( \phi\left( t \right) \right) \\ H_{phi(t)_{3}ˏ_3} =& \frac{\cos\left( \phi\left( t \right) \right)}{\cos\left( \theta\left( t \right) \right)} \\ \frac{\mathrm{d} \phi\left( t \right)}{\mathrm{d}t} =& H_{phi(t)_{1}ˏ_1} p\left( t \right) + H_{phi(t)_{1}ˏ_2} q\left( t \right) + H_{phi(t)_{1}ˏ_3} r\left( t \right) \\ \frac{\mathrm{d} \theta\left( t \right)}{\mathrm{d}t} =& H_{phi(t)_{2}ˏ_1} p\left( t \right) + H_{phi(t)_{2}ˏ_2} q\left( t \right) + H_{phi(t)_{2}ˏ_3} r\left( t \right) \\ \frac{\mathrm{d} \psi\left( t \right)}{\mathrm{d}t} =& H_{phi(t)_{3}ˏ_1} p\left( t \right) + H_{phi(t)_{3}ˏ_2} q\left( t \right) + H_{phi(t)_{3}ˏ_3} r\left( t \right) \end{align}

1.e. Optionally, a user may want to simulate the ODESystem. This step is not necessary for trimming the system but may provide intuition on the dynamics.

inputs = [uA, uT, uR, uE_1, uE_2]
outputs = []
sys, diff_idxs, alge_idxs, input_idxs = ModelingToolkit.io_preprocessing(rcam_model, inputs, outputs)

# Converting ODESystem to ODEProblem for numerical simulation.

x0 = Dict(
u => 85,
v => 0,
w => 0,
p => 0,
q => 0,
r => 0,
ϕ => 0,
θ => 0.1, # approx 5.73 degrees
ψ => 0
)

u0 = Dict(
uA => 0,
uT => -0.1,
uR => 0,
uE_1 => 0.08,
uE_2 => 0.08
)

tspan = (0.0, 3*60)
prob = ODEProblem(sys, x0, tspan, u0, jac = true)
sol = solve(prob, Tsit5())
plotvars = [u,
v,
w,
p,
q,
r,
ϕ,
θ,
ψ,]
1. Consider the case where the use desires to trimming the aircraft to steady-state straight-and-level flight. In the trim function, the default setting (when no derivative terms are specified in the arguments) automatically provides a steady-state solution. We assume that the user sets the airspeed Va to any arbitrary value (taken to be 85 [m/s] here). For straight-and-level flight, we assume that the flight path angle (the angle between the horizontal and velocity vector) given by the expression θ - atan(w,u) is set to 0 specified using a hard_eq_cons. We also set the following key-value pairs in the desired_states that v => 0 corresponding to no side-slip, ϕ => 0 to zero bank/roll angle and ψ => 0 to specify a North heading. We also specify inequality constraints using the hard_ineq_cons vector to denote bounds on the control inputs.
inputs = [uA, uT, uR, uE_1, uE_2]
outputs = []
trimmed_states = [u,
v,
w,
p,
q,
r,
ϕ,
θ,
ψ]
trimmed_controls = [uA,
uT,
uR,
uE_1,
uE_2]

desired_states = Dict(u => 85, v => 0, ϕ => 0, ψ => 0, uE_1 => 0.1, uE_2 => 0.1)

hard_eq_cons =  [
Va ~ 85
θ - atan(w,u) ~ 0.0
]

# These represent saturation limits on controllers.
hard_ineq_cons = [
]

penalty_multipliers = Dict(:desired_states => 10.0, :trim_cons => 10.0, :sys_eqns => 10.0)
sol, trim_states = trim(rcam_model; penalty_multipliers, desired_states, inputs, hard_eq_cons, hard_ineq_cons)
for elem in trimmed_states
println("$elem: ", round(trim_states[elem], digits = 4)) end for elem in trimmed_controls println("$elem:    ", round(trim_states[elem], digits = 4))
end
u(t):    84.9905
v(t):    -0.0
w(t):    1.2713
p(t):    0.0
q(t):    0.0
r(t):    -0.0
ϕ(t):    0.0
θ(t):    0.015
ψ(t):    0.0
uA(t):    0.0
uT(t):    -0.178
uR(t):    -0.0
uE_1(t):    0.0821
uE_2(t):    0.0821

The trim_states dictionary holds the values for the states and control inputs (including the auxiliary variables used to define the ODESystem prior to simplification). We may simulate the ODE using these initial values.

sys, diff_idxs, alge_idxs, input_idxs = ModelingToolkit.io_preprocessing(rcam_model, inputs, outputs)

# Converting ODESystem to ODEProblem for numerical simulation.

x0 = filter(trim_states) do (k,v)
k ∈ Set(trimmed_states)
end

u0 = filter(trim_states) do (k,v)
k ∈ Set(trimmed_controls)
end

tspan = (0.0, 3*60)
prob = ODEProblem(sys, x0, tspan, u0, jac = true)
sol = solve(prob, Tsit5())
plotvars = [u,
v,
w,
p,
q,
r,
ϕ,
θ,
ψ,]

plot(sol, idxs = plotvars, layout = length(plotvars), plot_title="Steady-state solution")

## DAE systems

Consider the semi-explicit DAE system described by the set of differential equations $\dot x = f(x,z,u,p)$ and algebraic equations $\mathbf{0} = g(x,z,u,p)$. The vector $x \in \mathbb{R}^n$ denotes the differential state variables (whose derivatives are present), $z \in \mathbb{R}^q$ denotes the algebraic variables (whose derivatives are not present), $u \in \mathbb{R}^m$ denotes control input variables, and $p$ denotes parameters. Functions $f:\mathbb{R}^{n} \times \mathbb{R}^{q} \times \mathbb{R}^{m} \mapsto \mathbb{R}^n$ and $g:\mathbb{R}^{n} \times \mathbb{R}^{q} \times \mathbb{R}^{m} \mapsto \mathbb{R}^q$ are potentially nonlinear functions. Similar to ODE systems, at an equilibrium point $(x_{eq}, z_{eq}, u_{eq})$ both $f(x_{eq}, z_{eq}, u_{eq}, p) = \mathbf{0}$ and $g(x_{eq}, z_{eq}, u_{eq}, p) = \mathbf{0}$ hold. At a trim point, the user may specify non-zero derivatives or impose additional desired constraints.

## Example: Electrical circuit with 2 inductors

Consider the circuit below where a voltage source is connected to a resistor and two inductors in series.

Figure 2: Electric circuit with voltage source, resistor and two inductors.

1.a. This is modeled with a component-wise approach using ModelingToolkitStandardLibrary. The resulting model has 28 equations and 29 variables.

using ModelingToolkit, JuliaSimControl, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks
using Plots

@parameters t

# Defining the parameters:
R = 4
L = 6
V = 12.0

# Defining the components and connections:
@named R1 = Resistor(; R = R)
@named L1 = Inductor(; L = 0.5*L)
@named L2 = Inductor(; L = 0.5*L )
@named voltage = Voltage()
@named input = RealInput()
@named output = RealOutput()
@named ground = Ground()
rl_source_eqns = [
Symbolics.connect(input, voltage.V)
Symbolics.connect(voltage.p, R1.p)
Symbolics.connect(R1.n, L1.p)
Symbolics.connect(L1.n, L2.p)
Symbolics.connect(L2.n, voltage.n)
Symbolics.connect(ground.g, voltage.n)
output.u ~ L2.i
]

# Defining the DAE system
@named circuit_dae = ODESystem(rl_source_eqns, t, systems=[R1, L1, L2, voltage, input, output, ground])
$$$$$\left[ \begin{array}{c} output_{+}u\left( t \right) = L2_{+}i\left( t \right) \\ \mathrm{connect}\left( input, voltage_{+}V \right) \\ \mathrm{connect}\left( voltage_{+}p, R1_{+}p \right) \\ \mathrm{connect}\left( R1_{+}n, L1_{+}p \right) \\ \mathrm{connect}\left( L1_{+}n, L2_{+}p \right) \\ \mathrm{connect}\left( L2_{+}n, voltage_{+}n \right) \\ \mathrm{connect}\left( ground_{+}g, voltage_{+}n \right) \\ R1_{+}v\left( t \right) = - R1_{+}n_{+}v\left( t \right) + R1_{+}p_{+}v\left( t \right) \\ 0 = R1_{+}n_{+}i\left( t \right) + R1_{+}p_{+}i\left( t \right) \\ R1_{+}i\left( t \right) = R1_{+}p_{+}i\left( t \right) \\ R1_{+}v\left( t \right) = R1_{+}R R1_{+}i\left( t \right) \\ L1_{+}v\left( t \right) = - L1_{+}n_{+}v\left( t \right) + L1_{+}p_{+}v\left( t \right) \\ 0 = L1_{+}n_{+}i\left( t \right) + L1_{+}p_{+}i\left( t \right) \\ L1_{+}i\left( t \right) = L1_{+}p_{+}i\left( t \right) \\ \frac{\mathrm{d} L1_{+}i\left( t \right)}{\mathrm{d}t} = \frac{L1_{+}v\left( t \right)}{L1_{+}L} \\ L2_{+}v\left( t \right) = - L2_{+}n_{+}v\left( t \right) + L2_{+}p_{+}v\left( t \right) \\ 0 = L2_{+}n_{+}i\left( t \right) + L2_{+}p_{+}i\left( t \right) \\ L2_{+}i\left( t \right) = L2_{+}p_{+}i\left( t \right) \\ \frac{\mathrm{d} L2_{+}i\left( t \right)}{\mathrm{d}t} = \frac{L2_{+}v\left( t \right)}{L2_{+}L} \\ voltage_{+}v\left( t \right) = - voltage_{+}n_{+}v\left( t \right) + voltage_{+}p_{+}v\left( t \right) \\ 0 = voltage_{+}n_{+}i\left( t \right) + voltage_{+}p_{+}i\left( t \right) \\ voltage_{+}i\left( t \right) = voltage_{+}p_{+}i\left( t \right) \\ voltage_{+}v\left( t \right) = voltage_{+}V_{+}u\left( t \right) \\ ground_{+}g_{+}v\left( t \right) = 0 \\ \end{array} \right]$$$$$

1.b. Optionally, the user may simplify the system. The voltage is set as a control input in the inputs array:

inputs = [input.u]
outputs = []
sys, diff_idxs, alge_idxs, input_idxs = ModelingToolkit.io_preprocessing(circuit_dae, inputs, outputs)
sys
\begin{align} \frac{\mathrm{d} L1_{+}i\left( t \right)}{\mathrm{d}t} =& L1_{+}iˍt\left( t \right) \\ 0 =& - L1_{+}iˍt\left( t \right) + \frac{L2_{+}v\left( t \right)}{L2_{+}L} \end{align}

After simplification, the system becomes a DAE with 1 algebraic equation, 1 ordinary differential equation, 2 states and 1 control input.

1.c. After simplification, one may simulate the DAE using a suitable algorithm such as Rodas4().

tspan = (0.0, 10)
x0 = Dict(
L1.p.i => 0.0,
L2.v => 0.0,
)
u0 = Dict(
input.u => 12
)
prob = ODEProblem(sys, x0, tspan, u0, jac = true)
sol = solve(prob, Rodas4())
plot(sol, idxs = [L1.i, L2.v])
1. To determine the trim point, here we set a desired current in the circuit, and call the trim function:
desired_states = Dict(L1.i => 9)
sol, trim_states = trim(circuit_dae; inputs, outputs, desired_states)
trim_states[L1.i]
8.999999989936743

## Index

JuliaSimControl.trimMethod
trim(sys::AbstractSystem;
inputs = [], outputs = [], params = Dict(),
desired_states = Dict(), hard_eq_cons = Equation[], soft_eq_cons = Equation[], hard_ineq_cons = Inequality[], soft_ineq_cons = Inequality[],
penalty_multipliers = Dict(:desired_states => 1.0, :trim_cons => 1.0, :sys_eqns => 10.0), dualize = false,
solver = IPNewton(), verbose = false, simplify = false, kwargs...)

Determine a trim point which is a feasible stationary point of the ODE closest to the specified desired state, given hard and soft constraints.

Arguments:

• sys::AbstractSystem: ODESystem to get a trim point for. Note that this should be the original system prior to callingstructural_simplify
• inputs = []: array of (control) inputs to the system
• outputs = []: array of outputs from the system.
• params: Dictionary of parameters that will override defaults.
• desired_states: Dictionary specifying the desired operating point or derivatives.
• hard_eq_cons = Equation[]: User-defined hard equality constraints.
• soft_eq_cons = Equation[]: User-defined soft equality constraints.
• hard_ineq_cons = Inequality[]: User-defined hard inequality constraints.
• soft_ineq_cons = Inequality[]: User-defined soft inequality constraints.
• penalty_multipliers::AbstractDict: The weighting of the :desired_states, :trim_cons (equality and inequality constraints), and :sys_eqns to the objective of the trimming optimization problem. Setting dualize to true and changing these multiples is useful in getting a trim point closer to user specifications.
• dualize = false: To dualize any user constraints as well as the system equations, set this to true. This solves an unconstrained optimization problem with the default solver changed to Newton().
• solver = IPNewton(): Specified optimization solver compatible with Optimization.jl. Note that the package for the interface to the solver needs to be loaded. By default only OptimizationOptimJL is added. In general, use a constrained optimization solver, so NOT Newton().
• verbose = false: Set to true for information prints.
• abstol = 0.0: User-defined absolute tolerance for constraints and ODESystem equations.
• simplify = false: Argument passed to ModelingToolkit.structural_simplify. Refer to documentation of ModelingToolkit.jl
• kwargs: Other keyword arguments passed to ModelingToolkit.find_solvables!

Note:

• If several hard and soft-constraints are imposed, it is important to give a good (feasible) initial point e.g., through desired for IPNewton() or set dualize to true.

Returns:

• 'sol::SciMLBase.OptimizationSolution': This provides information for the optimization solver. Refer to documentation of SciMLBase.jl
• 'trim_states: This provides a dictionary of all the states and derivatives ofsys as well as their corresponding values at the trim point.
source