Let \(s \in \mathbb{C}\). Use SymPy to perform a partial fraction expansion on the following expression: \[ \frac{(s+2) (s + 10)}{s^4 + 8 s^3 + 117 s^2 + 610 s + 500}. \]
First, load SymPy and define the expression:
import sympy as sp
s = sp.symbols("s", complex=True)
expr = (s + 2)*(s + 10)/(s**4 + 8*s**3 + 117*s**2 + 610*s + 500)The SymPy apart() method can
be directly applied to find the partial fraction expansion:
expr.apart()Let \(x, a_1, a_2, a_3, a_4 \in \mathbb{R}\). Use SymPy to combine the cosine and sine terms that share arguments into single sinusoids with phase shifts in the following expression: \[ a_1 \sin(x) + a_2 \cos(x) + a_3 \sin(2 x) + a_4 \cos(2 x) \]
First, load SymPy as follows:
import sympy as spUnfortunately, the SymPy trigsimp() method doesn’t combine the
cosine and sine terms that share an argument. However, we can define a
functions trig_two_to_one() for
performing the two-to-one conversion. The identities we will use are as
follows: $$\begin{align*}
a \cos u + b \sin u &= \sqrt{a^2 + b^2} \sin(u + \arctan(a/b)) \\
&= \sqrt{a^2 + b^2} \cos(u - \arctan(b/a)).
\end{align*}$$ We begin by defining two functions, each of which returns
a dictionary of replacement rules for sin(u) in the identities
above.
def _trig_two_to_one_sin_rule(a, b, u):
"""Replacement rule for sin(u) in a cos(u) + b sin(u) for
converting to a single sin term with a phase
The identity applied is:
a cos u + b sin u = sqrt(a^2 + b^2) sin(u + atan(a/b))
Returns: A dictionary for replacing sin(u)
"""
identity = a*sp.cos(u) + b*sp.sin(u) + \
- sp.sqrt(a**2 + b**2)*sp.sin(u + sp.atan(a/b))
sol = sp.solve(identity, [sp.sin(u)], dict=True)
return sol[0]
def _trig_two_to_one_cos_rule(a, b, u):
"""Replacement rule for sin(u) in a cos(u) + b sin(u) for
converting to a single cos term with a phase
The identity applied is:
a cos u + b sin u = sqrt(a^2 + b^2) cos(u - atan(b/a))
Returns: A dictionary for replacing sin(u)
"""
identity = a*sp.cos(u) + b*sp.sin(u) + \
- sp.sqrt(a**2 + b**2)*sp.cos(u - sp.atan(b/a))
sol = sp.solve(identity, [sp.sin(u)], dict=True)
return sol[0]Now we can write a function to perform the trigonometric simplification, as follows:
def trig_two_to_one(expr: sp.Expr, to: str = "sin"):
"""Rewrites sin and cos terms that share arguments to single sin
or cos terms
Applies the following identity:
a cos u + b sin u = sqrt(a^2 + b^2) sin(u + atan(a/b))
or
a cos u + b sin u = sqrt(a^2 + b^2) cos(u - atan(b/a))
depending on the ``to`` argument.
Args:
expr: The symbolic expression containing sin(u) and cos(u)
to: Rewrite with "sin" (default) or "cos"
"""
expr = expr.simplify()
# Identify sin terms
w1 = sp.Wild("w1", exclude=[1])
w2 = sp.Wild("w2")
sin_terms = expr.find(w1*sp.sin(w2))
sin_arg_amps = {} # To be: {sin argument: sine amplitude}
for term in sin_terms:
arg_amp_rules = term.match(w1*sp.sin(w2))
sin_arg_amps[arg_amp_rules[w2]] = arg_amp_rules[w1]
# Identify cos terms
cos_terms = expr.find(w1*sp.cos(w2))
cos_arg_amps = {} # To be: {sin argument: sine amplitude}
for term in cos_terms:
arg_amp_rules = term.match(w1*sp.cos(w2))
cos_arg_amps[arg_amp_rules[w2]] = arg_amp_rules[w1]
# Replace with wildcard rule
for sin_arg, sin_amp in sin_arg_amps.items():
if to == "sin":
if sin_arg in cos_arg_amps.keys():
cos_amp = cos_arg_amps[sin_arg]
sin_rule = _trig_two_to_one_sin_rule(
cos_amp, sin_amp, sin_arg
)
expr = expr.subs(sin_rule)
elif to == "cos":
if sin_arg in cos_arg_amps.keys():
cos_amp = cos_arg_amps[sin_arg]
cos_rule = _trig_two_to_one_cos_rule(
cos_amp, sin_amp, sin_arg
)
expr = expr.subs(cos_rule)
return expr.simplify()The expression of interest in the problem is defined here:
x, a1, a2, a3, a4 = sp.symbols("x, a1, a2, a3, a4", real=True)
expr = a1*sp.sin(x) + a2*sp.cos(x) + a3*sp.sin(2*x) + a4*sp.cos(2*x)Apply the function to the expression, rewriting in terms of sine (first) and cosine (second):
trig_two_to_one(expr, to="sin")
trig_two_to_one(expr, to="cos")Consider the following equation, where \(x \in \mathbb{C}\) and \(a, b, c \in \mathbb{R_+}\), \[ a x^2 + b x + \frac{c}{x} + b^2 = 0. \] Use SymPy to solve for \(x\).
First, load SymPy as follows:
import sympy as spNow define the symbolic variables and the equation:
x = sp.symbols("x", complex=True)
a, b, c = sp.symbols("a, b, c", positive=True)
eq = a*x**2 + b*x + c/x + b**2 # == 0Let’s try sp.solve(), as
follows:
sol = sp.solve(eq, x, dict=True)
print(sol)[{x: -(-3*b**2/a + b**2/a**2)/(3*(sqrt(-4*(-3*b**2/a + b**2/a**2)**3 + (27*c/a - 9*b**3/a**2 + 2*b**3/a**3)**2)/2 + 27*c/(2*a) - 9*b**3/(2*a**2) + b**3/a**3)**(1/3)) - (sqrt(-4*(-3*b**2/a + b**2/a**2)**3 + (27*c/a - 9*b**3/a**2 + 2*b**3/a**3)**2)/2 + 27*c/(2*a) - 9*b**3/(2*a**2) + b**3/a**3)**(1/3)/3 - b/(3*a)}, {x: -(-3*b**2/a + b**2/a**2)/(3*(-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*b**2/a + b**2/a**2)**3 + (27*c/a - 9*b**3/a**2 + 2*b**3/a**3)**2)/2 + 27*c/(2*a) - 9*b**3/(2*a**2) + b**3/a**3)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*b**2/a + b**2/a**2)**3 + (27*c/a - 9*b**3/a**2 + 2*b**3/a**3)**2)/2 + 27*c/(2*a) - 9*b**3/(2*a**2) + b**3/a**3)**(1/3)/3 - b/(3*a)}, {x: -(-3*b**2/a + b**2/a**2)/(3*(-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*b**2/a + b**2/a**2)**3 + (27*c/a - 9*b**3/a**2 + 2*b**3/a**3)**2)/2 + 27*c/(2*a) - 9*b**3/(2*a**2) + b**3/a**3)**(1/3)) - (-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*b**2/a + b**2/a**2)**3 + (27*c/a - 9*b**3/a**2 + 2*b**3/a**3)**2)/2 + 27*c/(2*a) - 9*b**3/(2*a**2) + b**3/a**3)**(1/3)/3 - b/(3*a)}]This is a cubic solution, so it is unwieldy.
Let \(w, x, y, z \in \mathbb{R}\). Consider the following system of equations: $$\begin{alignat*}{3} 8 w &{}-{} 6 x &{}+{} 5 y &{}+{} 4 z &{}={} -20 \\ \phantom{3 w} &{}\mathbin{\phantom{+}}{} \phantom{4 x} &{}\mathbin{\phantom{+}}{} 2 y &{}-{} 2 z &{}={} 10 \\ 2 w &{}-{} \phantom{3}x &{}+{} 4 y &{}+{} \phantom{3}z &{}={} 0\\ w &{}+{} 4 x &{}-{} 2 y &{}+{} 8 z &{}={} 4. \end{alignat*}$$ Use SymPy to solve the system for \(w\), \(x\), \(y\), and \(z\).
Load SymPy as follows:
import sympy as spNow define the symbolic variables:
w, x, y, z = sp.symbols("w, x, y, z", complex=True)Define the set of equations:
S = [
8*w - 6*x + 5*y + 4*z + 20, # == 0
2*y - 2*z - 10, # == 0
2*w - x + 4*y + z, # == 0
w + 4*x - 2*y + 8*z - 4, # == 0
]Let’s try sp.solve(), as
follows:
sol = sp.solve(S, [w, x, y, z], dict=True)
print(sol)[{w: 564/85, x: 773/85, y: 14/85, z: -411/85}]Consider the truss shown in fig. ¿fig:truss-03?. Use a static analysis and the method of joints to develop a solution for the force in each member \(F_{AC}\), \(F_{AD}\), etc., and the reaction forces using the sign convention that tension is positive and compression is negative. The forces should be expressed in terms of the applied force \(\vec{f}_D\) and the dimensions \(w\) and \(h\) only. Write a program that solves for the forces symbolically and answers the following questions:
- Which members are in tension?
- Which members are in compression?
- Are there any members with \(0\) nominal force? If so, which?
- Which member (or members) has (or have) the maximum compression?
- Which member (or members) has (or have) the maximum tension?
#### Import Packages {-}import numpy as np
import sympy as sp
import matplotlib.pyplot as pltStatic Analysis
Using the method of joints, we proceed through the joints, summing forces in the \(x\)- and \(y\)-directions. We will assume all members are in tension, and their sign will be positive if this is the case and negative, otherwise. Beginning with joint A, which includes two reaction forces \(R_{A_x}\) and \(R_{A_y}\) from the support, $$\begin{align} \Sigma F_x &= 0; & R_{A_x} + F_{AD} + F_{AC} \cos\theta &=0 \\ \Sigma F_y &= 0; & R_{A_y} - F_{AC} \sin\theta &= 0. \end{align}$$ The angle \(\theta\) is known in terms of the dimensions \(w\) and \(h\) as \[\theta = \arctan\frac{h}{w}.\] These equations can be encode symbolically as follows:
RAx, RAy, FAC, FAD, theta = sp.symbols(
"RAx, RAy, FAC, FAD, theta", real=True
)
h, w = sp.symbols("h, w", positive=True)
eqAx = RAx + FAD + FAC * sp.cos(theta)
eqAy = RAy - FAC*sp.sin(theta)
theta_wh = sp.atan(h/w)Proceeding to joint B, which also has two reaction forces \(R_{B_x}\) and \(R_{B_y}\), $$\begin{align} \Sigma F_x &= 0; & R_{B_x} + F_{BC} &= 0 \\ \Sigma F_y &= 0; & R_{B_y} &= 0. \end{align}$$ Encoding these equations,
RBx, RBy, FBC = sp.symbols("RBx, RBy, FBC", real=True)
eqBx = RBx + FBC
eqBy = RByFor joint C, $$\begin{align} \Sigma F_x &= 0; & -F_{BC} - F_{AC} \cos\theta &= 0 \\ \Sigma F_y &= 0; & F_{CD} + F_{AC} \sin\theta &= 0. \end{align}$$ Encoding these equations,
FCD = sp.symbols("FCD", real=True)
eqCx = -FBC - FAC * sp.cos(theta)
eqCy = FCD + FAC * sp.sin(theta)For joint D, there is an applied load \(\vec{f_D}\), so the analysis proceeds as follows: $$\begin{align} \Sigma F_x &= 0; & -F_{AD} - f_D \cos\frac{\pi}{4} &= 0 \\ \Sigma F_y &= 0; & -F_{CD} - f_D \sin\frac{\pi}{4} &= 0. \end{align}$$ Encoding these equations,
fD = sp.symbols("fD", positive=True)
eqDx = -FAD - fD * sp.cos(sp.pi/4)
eqDy = -FCD - fD * sp.sin(sp.pi/4)In total, we have 8 force equations and 8 unknown forces (4 member forces and 4 reaction forces). Let’s construct the system and solve it for the unkown forces, as follows:
S_forces = [
eqAx, eqAy, eqBx, eqBy, eqCx, eqCy, eqDx, eqDy
] # 8 force equations
member_forces = [FAC, FAD, FBC, FCD] # 5 member forces
reaction_forces = [RAx, RAy, RBx, RBy] # 3 reaction forces
forces_unknown = member_forces + reaction_forces # 8 unkown forces
sol_forces = sp.solve(S_forces, forces_unknown, dict=True); sol_forces[{FAC: sqrt(2)*fD/(2*sin(theta)),
FAD: -sqrt(2)*fD/2,
FBC: -sqrt(2)*fD*cos(theta)/(2*sin(theta)),
FCD: -sqrt(2)*fD/2,
RAx: sqrt(2)*fD/2 - sqrt(2)*fD*cos(theta)/(2*sin(theta)),
RAy: sqrt(2)*fD/2,
RBx: sqrt(2)*fD*cos(theta)/(2*sin(theta)),
RBy: 0}]This solution is in terms of \(f_D\)
and \(\theta\). Because \(w\) and \(h\) are our design parameters, let’s
substitute theta_wh such that our
solution is rewritten in terms of \(f_D\), \(w\), and \(h\). Create a dict of
solutions as follows:
forces_wh = {} # Initialize
for force in forces_unknown:
force_wh = force.subs(
sol_forces[0]
).subs(
theta, theta_wh
).simplify()
forces_wh[force] = force_wh
print(f"{force} = {force_wh}")FAC = fD*sqrt(2*h**2 + 2*w**2)/(2*h)
FAD = -sqrt(2)*fD/2
FBC = -sqrt(2)*fD*w/(2*h)
FCD = -sqrt(2)*fD/2
RAx = sqrt(2)*fD*(h - w)/(2*h)
RAy = sqrt(2)*fD/2
RBx = sqrt(2)*fD*w/(2*h)
RBy = 0Identify the member forces in tension and compression as follows:
tension = [] # Initialize list
compression = [] # Initialize list
for force in member_forces:
if forces_wh[force] > 0:
tension.append(force)
print(f"{force} is in tension")
elif forces_wh[force] < 0:
compression.append(force)
print(f"{force} is in compression")
else:
print(f"{force} is a zero-force member")FAC is in tension
FAD is in compression
FBC is in compression
FCD is in compressionTherefore, there are no zero-force members. The maximum compression depends on the ratio \(w/h\). If \(w/h = 1\), all three members in compression have the same compression of magnitude \[ \frac{\sqrt{2}}{2} f_D. \] This is also the maximum compression if \(w/h < 1\). If \(w/h > 1\), the maximum compression is of magnitude \[ \frac{\sqrt{2} w}{2 h} f_D. \]
You are designing the truss structure shown in fig. ¿fig:truss-02?, which is to support the hanging of an external load \(\bm{f}_C = -f_C \bm{\hat{j}}\), where \(f_C > 0\). Your organization plans to offer customers the following options:
- Any width (i.e., \(2 w\))
- A selection of maximum load magnitudes \(L = f_C / \alpha \in \Gamma\), where \(\Gamma = \{1, 2, 4, 8, 16\}\) kN, and where \(\alpha\) is the factor of safety
As the designer, you are to develop a design curve for the dimension \(h\) versus half-width \(w\) for each maximum load \(L \in \Gamma\), under the following design constraints:
- Minimize the dimension \(h\)
- The tension in all members is no more than a given \(T\)
- The compression in all members is no more than a given \(C\)
- The magnitude of the support force at pin A is no more than a given \(P_A\)
- The magnitude of the support force at pin D is no more than a given \(P_D\)
Use a static analysis and the method of joints to develop a solution for the force in each member \(F_{AB}\), \(F_{AC}\), etc., and the reaction forces using the sign convention that tension is positive and compression is negative. Create a Python function that returns \(h\) as a function of \(w\) for a given set of design parameters \(\{T, C, P_A, P_D, \alpha, L\}\). Use the function to create a design curve \(h\) versus \(2 w\) for each \(L \in \Gamma\), maximum tension \(T = 81\) kN, maximum compression \(C = 81\) kN, maximum support A load \(P_A = 50\) kN, maximum support D load \(P_D = 50\) kN, and a factor of safety of \(\alpha = 5\).
#### Import Packages {-}import numpy as np
import sympy as sp
import matplotlib.pyplot as pltStatic Analysis
Using the method of joints, we proceed through the joints, summing forces in the \(x\)- and \(y\)-directions. We will assume all members are in tension, and their sign will be positive if this is the case and negative, otherwise. Beginning with joint A, which includes one reaction force \(R_{A}\) from the support, $$\begin{align} \Sigma F_x &= 0; & F_{AC} + F_{AB} \cos\theta &=0 \\ \Sigma F_y &= 0; & R_A + F_{AB} \sin\theta &= 0. \end{align}$$ The angle \(\theta\) is known in terms of the dimensions \(w\) and \(h\) as \[\theta = \arctan\frac{h}{w}.\] These equations can be encode symbolically as follows:
RA, FAB, FAC, theta= sp.symbols(
"RA, FAB, FAC, theta", real=True
)
h, w = sp.symbols("h, w", positive=True)
eqAx = FAC + FAB*sp.cos(theta)
eqAy = RA + FAB*sp.sin(theta)
theta_wh = sp.atan(h/w)Proceeding to joint B, $$\begin{align} \Sigma F_x &= 0; & -F_{AB} \cos\theta + F_{BD} \cos\theta &= 0 \\ \Sigma F_y &= 0; & -F_{AB} \sin\theta - F_{BD} \sin\theta - F_{BC} &= 0. \end{align}$$ Encoding these equations,
FBD, FBC = sp.symbols("FBD, FBC", real=True)
eqBx = -FAB*sp.cos(theta) + FBD*sp.cos(theta)
eqBy = -FAB*sp.sin(theta) - FBD*sp.sin(theta) - FBCFor joint C, there is an externally applied force \(\bm{f}_C\), so the analysis proceeds as follows: $$\begin{align} \Sigma F_x &= 0; & -F_{AC} + F_{CD} &= 0 \\ \Sigma F_y &= 0; & F_{BC} - f_C &= 0. \end{align}$$ Encoding these equations,
FCD = sp.symbols("FCD", real=True)
fC = sp.symbols("fC", positive=True)
eqCx = -FAC + FCD
eqCy = FBC - fCFor joint D, the pinned reaction forces \(R_{Dx}\) and \(R_{Dy}\) are present, so the analysis proceeds as follows: $$\begin{align} \Sigma F_x &= 0; & -F_{CD} - F_{BD} \cos\theta + R_{Dx} &= 0 \\ \Sigma F_y &= 0; & F_{BD} \sin\theta + R_{Dy} &= 0. \end{align}$$ Encoding these equations,
RDx, RDy = sp.symbols("RDx, RDy", real=True)
eqDx = -FCD - FBD*sp.cos(theta) + RDx
eqDy = FBD*sp.sin(theta) + RDyIn total, we have 8 force equations and 8 unknown forces (5 member forces and 3 reaction forces). Let’s construct the system and solve it for the unkown forces, as follows:
S_forces = [
eqAx, eqAy, eqBx, eqBy, eqCx, eqCy, eqDx, eqDy
] # 8 force equations
member_forces = [FAB, FAC, FBC, FBD, FCD] # 5 member forces
reaction_forces = [RA, RDx, RDy] # 3 reaction forces
forces_unknown = member_forces + reaction_forces # 8 unkown forces
sol_forces = sp.solve(S_forces, forces_unknown, dict=True); sol_forces[{FAB: -fC/(2*sin(theta)),
FAC: fC*cos(theta)/(2*sin(theta)),
FBC: fC,
FBD: -fC/(2*sin(theta)),
FCD: fC*cos(theta)/(2*sin(theta)),
RA: fC/2,
RDx: 0,
RDy: fC/2}]This solution is in terms of \(f_C\)
and \(\theta\). Because \(w\) and \(h\) are our design parameters, let’s
substitute eqtheta such that our
solution is rewritten in terms of \(f_F\), \(w\), and \(h\). Create a list of solutions as
follows:
forces_wh = {} # Initialize
for force in forces_unknown:
force_wh = force.subs(
sol_forces[0]
).subs(
theta, theta_wh
).simplify()
forces_wh[force] = force_wh
print(f"{force} = {force_wh}")FAB = -fC*sqrt(h**2 + w**2)/(2*h)
FAC = fC*w/(2*h)
FBC = fC
FBD = -fC*sqrt(h**2 + w**2)/(2*h)
FCD = fC*w/(2*h)
RA = fC/2
RDx = 0
RDy = fC/2Define the constraints as follows:
C, T, PA, PD = sp.symbols("C, T, PA, PD", positive=True)
member_constraints = {"Tension": T, "Compression": C}
reaction_constraints = {
PA: (RA, 0),
PD: (RDx, RDy)
} # Max magnitude: (x force, y force)Define a function that encodes the constraints as a list of expressions that must be nonnegative.
def get_force_constraints(
forces: dict,
member_forces: list,
member_constraints: dict,
reaction_constraints: dict,
):
"""Returns a list of expressions that must be nonnegative
Args:
- forces: Force solutions {force symbol: solution}
- member_forces: List of member symbols [FAB, FAC ...]
- member_constraints:
{"Tension": max tension, "Compression": max compression}
- reaction_constraints:
{max force symbol: (max x force, max y force)}
"""
force_constraints = []
# Append member constraints
for f_name, f_value in forces.items():
if f_name in member_forces:
if f_value > 0: # Tension
force_constraints.append(
member_constraints["Tension"] - f_value
)
elif f_value < 0: # Compression
force_constraints.append(
member_constraints["Compression"] + f_value
)
# Append reaction constraints
for constraint, pair in reaction_constraints.items():
force_constraints.append(
constraint - sp.sqrt(pair[0]**2 + pair[1]**2).subs(forces_wh)
)
return force_constraintsNow apply the function:
constraints = get_force_constraints(
forces_wh, member_forces, member_constraints, reaction_constraints
)
print(constraints)[C - fC*sqrt(h**2 + w**2)/(2*h), T - fC*w/(2*h), T - fC, C - fC*sqrt(h**2 + w**2)/(2*h), T - fC*w/(2*h), PA - fC/2, PD - fC/2]At this point, we can solve for \(h\) in each expression that includes \(h\). Because these constraints are to be nonnegative, solving with equality gives the minimum \(h_\text{min}\). Expressions that don’t include \(h\) must still be checked. We can use this opportunity to sort the constraint expressions into those that involve \(h\) and those that do not.
h_mins = []
constraints_to_check = []
for constraint in constraints:
if h in constraint.free_symbols:
# Solve for h
h_min_sol = sp.solve(constraint, h)
for h_min in h_min_sol:
h_mins.append(h_min)
else:
constraints_to_check.append(constraint)Inspect the h_mins and constraints_to_check:
print(f"Min h solutions: {h_mins}")
print(f"Constraints to be checked: {constraints_to_check}")Min h solutions: [fC*w*sqrt(1/(2*C - fC))/sqrt(2*C + fC), fC*w/(2*T), fC*w*sqrt(1/(2*C - fC))/sqrt(2*C + fC), fC*w/(2*T)]
Constraints to be checked: [T - fC, PA - fC/2, PD - fC/2]Finally, define a function to design the truss:
def truss_designer(
h_mins: list,
constraints_to_check: list,
design_params: dict,
):
"""Returns an expression for h(w) using the design parameters provided
"""
L = sp.symbols("L", positive=True)
alpha = sp.symbols("alpha", positive=True)
for constraint in constraints_to_check:
constraint_ = constraint.subs(fC, L*alpha).subs(design_params)
if constraint_ < 0:
raise Exception(
f"Design failed for constraint: {constraint} < 0"
)
h_mins_ = []
for h_min in h_mins:
h_mins_.append(h_min.subs(fC, L*alpha).subs(design_params))
return max(h_mins_) # Maximum h_minDefine the given design parameters in a dictionary:
alpha = sp.symbols("alpha", positive=True)
design_parameters = {
T: 81e3, C: 81e3, PA: 50e3, PD: 50e3, alpha: 5,
} # Forces in N
Gamma = [1e3, 2e3, 4e3, 8e3, 16e3]Loop through the L loads, running
truss_designer():
h_mins_w = []
L = sp.symbols("L", positive=True)
for L_ in Gamma:
design_parameters[L] = L_
h_min = truss_designer(
h_mins, constraints_to_check, design_parameters
)
h_mins_w.append(h_min)
print(h_min)0.0308789086390917*w
0.0618463369994117*w
0.124408521678431*w
0.254802914503365*w
0.56790458868584*wConvert these symbolic expressions into numerically evaluable function, as follows:
h_min_funs = []
for h_min in h_mins_w:
h_min_funs.append(
sp.lambdify([w], h_min, modules=np)
)Plot the design curves as follows:
w_ = np.linspace(1, 20, 101)
fig, ax, = plt.subplots()
for ih, h_min_fun in enumerate(h_min_funs):
ax.semilogy(2*w_, h_min_fun(w_), label=f"L = {Gamma[ih]}")
ax.set_xlabel(f"width $2w$ (m)")
ax.set_ylabel(f"height $h$ (m)")
ax.grid()
ax.legend()
plt.show()
Consider an LTI system modeled by the state equation of the state-space model, eq. ¿eq:state-space-model-state?. A steady state of a system is defined as the state vector \(\bm{x}(t)\) after the effects of initial conditions have become relatively small. For a constant input \(\bm{u}(t) = \overline{\bm{u}}\), the constant state \(\overline{\bm{x}}\) toward which the system’s response decays can be found by setting the time derivative vector \(\bm{x}'(t) = \bm{0}\).
Write a Python function steady_state() that accepts the
following arguments:
A: A symbolic matrix representing \(A\)B: A symbolic matrix representing \(B\)u_const: A symbolic vector representing \(\overline{\bm{u}}\)
The function should return x_const, a symbolic vector representing
\(\overline{\bm{x}}\).
The steady-state output converges to \(\overline{\bm{y}}\) the corresponding
output equation of the state-space model,
eq. ¿eq:state-space-model-output?. Write a second
Python function steady_output()
that accepts the following arguments:
C: A symbolic matrix representing \(C\)D: A symbolic matrix representing \(D\)u_const: A symbolic vector representing \(\overline{\bm{u}}\)x_const: A symbolic vector representing \(\overline{\bm{x}}\)
This function should return y_const, a symbolic vector representing
\(\overline{\bm{y}}\).
Apply steady_state() and steady_output() to the state-space
model of the circuit shown in fig. ¿fig:w5?, which
includes a resistor with resistance \(R\), an inductor with inductance \(L\), and capacitor with capacitance \(C\). The LTI system is represented by
eq. ¿eq:state-space-model? with state, input, and
output vectors
\[
\bm{x}(t) = \begin{bmatrix} v_C(t) \\ i_L(t) \end{bmatrix},\
\bm{u}(t) = \begin{bmatrix} V_S \end{bmatrix},\
\bm{y}(t) = \begin{bmatrix} v_C(t) \\ v_L(t) \end{bmatrix}
\] and the following matrices: \[
A = \begin{bmatrix} 0 & 1/C \\ -1/L & -R/L \end{bmatrix},\
B = \begin{bmatrix} 0 \\ 1/L \end{bmatrix},\
C = \begin{bmatrix} 1 & 0 \\ -1 & -R \end{bmatrix},\
D = \begin{bmatrix} 0 \\ 1 \end{bmatrix}.
\] Furthermore, let the constant input vector be \[
\overline{\bm{u}} = \begin{bmatrix} \overline{V_S} \end{bmatrix},
\] for constant \(\overline{V_S}\).
A constant steady-state, \(\bm{x}' =
\bm{0}\) implies, from the state
eq. ¿eq:state-space-model-state?, $$\begin{align}
\bm{0} &= A \overline{\bm{x}} + B \overline{\bm{u}} \Rightarrow \\
\overline{\bm{x}} &= -A^{-1} B \overline{\bm{u}}.
\end{align}$$ We are now ready to define steady_state() as follows:
def steady_state(A, B, u_const):
"""Returns the symbolic constant steady state vector"""
A = sp.Matrix(A) # In case A isn't symbolic
B = sp.Matrix(B) # In case B isn't symbolic
u_const = sp.Matrix(u_const) # In case u_const isn't symbolic
x_const = -A**-1 * B * u_const
return x_constThe state-space output equation
eq. ¿eq:state-space-model-output? is already solved for
the output, so we are ready to write steady_output() as follows:
def steady_output(C, D, u_const, x_const):
"""Returns the symbolic constant steady-state output vector"""
C = sp.Matrix(C) # In case C isn't symbolic
D = sp.Matrix(D) # In case D isn't symbolic
u_const = sp.Matrix(u_const) # In case u_const isn't symbolic
x_const = sp.Matrix(x_const) # In case x_const isn't symbolic
y_const = C*x_const + D*u_const
return y_constApply these functions to the given state-space model. First, define the symbolic variables as follows:
R, L, C1 = sp.symbols("R, L, C1", positive=True)
VS_ = sp.symbols("VS_", real=True) # Constant voltage source inputNow define the system and the constant input as follows:
A = sp.Matrix([[0, 1/C1], [-1/L, -R/L]]) # $A$
B = sp.Matrix([[0], [1/L]]) # $B$
C = sp.Matrix([[1, 0], [-1, -R]]) # $C$
D = sp.Matrix([[0], [1]]) # $D$
u_const = sp.Matrix([[VS_]]) # $\overline{\bm{u}}$Find the constant steady state \(\overline{\bm{x}}\) as follows:
x_const = steady_state(A, B, u_const)
print(x_const)Find the constant steady-state output \(\overline{\bm{y}}\) as follows:
y_const = steady_output(C, D, u_const, x_const)
print(y_const)Consider the electromechanical state-space model described in example. For a given set of parameters, input voltage, and initial conditions, the following vector-valued functions have been derived: \[ \bm{F} = \begin{bmatrix} \int_0^t v_R(t)\ dt \\ \int_0^t v_L(t)\ dt \\ \int_0^t \Omega_B(t)\ dt \\ \int_0^t \Omega_J(t)\ dt \\ \end{bmatrix} = \begin{bmatrix} \exp(-t) \\ \exp(-t) \\ 1 - \exp(-t) \\ 1 - \exp(-t) \end{bmatrix}, \quad \bm{G} = \begin{bmatrix} \int_0^t i_R(t)\ dt \\ \int_0^t i_L(t)\ dt \\ \int_0^t T_B(t)\ dt \\ \int_0^t T_J(t)\ dt \\ \end{bmatrix} = \begin{bmatrix} \exp(-t) \\ \exp(-t) \\ 1 - \exp(-t) \\ \exp(-t) \end{bmatrix} \] The instantaneous power lossed or stored by each element is given by the following vector of products: \[ \bm{\mathcal{P}}(t) = \begin{bmatrix} v_R(t) i_R(t) \\ v_L(t) i_L(t) \\ \Omega_B(t) T_B(t) \\ \Omega_J(t) T_J(t) \end{bmatrix}. \] The energy \(\mathcal{E}(t)\) of the elements, then, is \[ \bm{\mathcal{E}}(t) = \int_0^t \bm{\mathcal{P}}(t) dt. \]
Write a program that satisfies the following requirements:
- It defines a function
power(F, G)that returns the symbolic power vector \(\bm{\mathcal{P}}(t)\) from any inputs \(\bm{F}\) and \(\bm{G}\) - It defines a function
energy(F, G)that returns the symbolic energy \(\bm{\mathcal{E}}(t)\) from any inputs \(\bm{F}\) and \(\bm{G}\) (energy()should callpower()) - It tests the
energy()on the specific \(\bm{F}\) and \(\bm{G}\) given above
The formula for the power of each element is given, so we are ready
to define power() as follows:
def power(F, G):
"""Returns the power for vectors F and G"""
F = sp.Matrix(F) # In case F isn't symbolic
G = sp.Matrix(G) # In case G isn't symbolic
P = F.multiply_elementwise(G)
# Alternative using a for loop:
# P = sp.zeros(*F.shape) # Initialize
# for i, Fi in enumerate(F):
# P[i] = Fi * G[i]
return PThe formula for the energy stored or dissipated by each element is
given, so we are ready to write energy() as follows:
def energy(F, G):
"""Returns the energy stored for vectors F and G"""
P = power(F, G)
E = sp.integrate(P, (t, 0, t))
return EApply these functions to the given \(\bm{F}\) and \(\bm{G}\). First, define \(\bm{F}\) and \(\bm{G}\) as follows:
t = sp.symbols("t", real=True)
F = sp.Matrix([
[sp.exp(-t)],
[sp.exp(-t)],
[1 - sp.exp(-t)],
[1 - sp.exp(-t)]
])
G = sp.Matrix([
[sp.exp(-t)],
[sp.exp(-t)],
[1 - sp.exp(-t)],
[sp.exp(-t)]
])Now compute the energy:
E = energy(F, G).simplify()
print(E)For the circuit and state-space model given in problem 4.7, use SymPy to solve for \(\vec{x}(t)\) and \(\vec{y}(t)\) given the following:
- A constant input voltage \(V_S(t) = \overline{V_S}\)
- Initial condition \(\vec{x}(0) = \vec{0}\)
Substitute the following parameters into the solution for \(\vec{y}(t)\) and create numerically evaluable functions of time for each variable in \(\vec{y}(t)\): \[ R = 50\text{ $\Omega$, } L = 10\cdot 10^{-6}\text{ H, } C = 1\cdot 10^{-9}\text{ F, } \overline{V_S} = 10\text{ V}. \]
Plot the outputs in \(\vec{y}(t)\) as functions of time, making sure to choose a range of time over which the response is best presented. Hint: An appropriate amount of time is on the scale of microseconds.
Load packages:
import numpy as np
import sympy as sp
import matplotlib.pyplot as pltDefine Classes
We begin by defining the parameters and functions of time as SymPy symbolic variables and unspecified functions as follows:
R, L, C = sp.symbols("R, L, C", positive=True)
v_C, i_L, v_L, V_S = sp.symbols(
"v_C, i_L, v_L, V_S", cls=sp.Function, real=True
) # $v_C, i_L, V_S$
t = sp.symbols("t", real=True)Now we can form the symbolic matrices and vectors:
A_ = sp.Matrix([[0, 1/C], [-1/L, -R/L]]) # $\mat{A}$
B_ = sp.Matrix([[0], [1/L]]) # $\mat{B}$
C_ = sp.Matrix([[1, 0], [-1, -R]]) # $\mat{C}$
D_ = sp.Matrix([[0], [1]]) # $\mat{D}$
x = sp.Matrix([[v_C(t)], [i_L(t)]]) # $\vec{x}$
u = sp.Matrix([[V_S(t)]]) # $\vec{u}$
y = sp.Matrix([[v_C(t)], [v_L(t)]]) # $\vec{y}$The input and initial conditions can be encoded as follows:
u_subs = {V_S(t): 10}
ics = {v_C(0): 0, i_L(0): 0}The set of first-order ODEs comprising the state equation can be defined as follows:
odes = x.diff(t) - A_*x - B_*u
print(odes)x_sol = sp.dsolve(list(odes.subs(u_subs)), list(x), ics=ics)The symbolic solutions for \(\bm{x}(t)\) are lengthy expressions, so we don’t print them here. Now we can compute the output \(\vec{y}(t)\) from eq. ¿eq:state-space-model-output? as follows:
x_sol_dict = {} # Initialize
for eq in x_sol:
x_sol_dict[eq.lhs] = eq.rhs # Make a dict of solutions for subs
y_sol = (C_*x + D_*u).subs(x_sol_dict) # Subs into output equation
# We will graph the output for the following set of parameters:params = {
R: 50, # (Ohms)
L: 10e-6, # (H)
C: 1e-9, # (F)
}Create a numerically evaluable version of each function as follows:
v_C_ = sp.lambdify(
t, y_sol[0].subs(params).subs(u_subs), modules="numpy"
)
v_L_ = sp.lambdify(
t, y_sol[1].subs(params).subs(u_subs), modules="numpy"
)Graph each solution as follows:
t_ = np.linspace(0, 0.000002, 201)
fig, axs = plt.subplots(2, sharex=True)
axs[0].plot(t_, v_C_(t_))
axs[1].plot(t_, v_L_(t_))
axs[1].set_xlabel("Time (s)")
axs[0].set_ylabel("$v_C(t)$ (rad/s)")
axs[1].set_ylabel("$v_L(t)$ (A)")
plt.show()
The output equation is trivial in this case, yielding only the state variable \(\Omega_J(t)\), for which we have already solved. Therefore, we have completed the analysis.
Online Resources for Section 4.8
No online resources.