In [6]:
#question 1
import numpy as np
import [Link] as plt
import control as ctrl
G = [Link]([25], [1, 5, 0])
T = [Link](G, 1, -1)
print("Part a - Closed-loop TF:\n", T)
poles = [Link](T)
print("Part b - Poles:", poles)
t, y = ctrl.step_response(T)
[Link](t, y)
[Link]("Part c: Step Response")
[Link]()
[Link]()
print("performance metrics", ctrl.step_info(T))
mag1, ph1, w1 = ctrl.frequency_response(G, 1)
print("mag at 1 rad/s:", mag1, "Phase:", ph1)
mag10, ph10, w10 = ctrl.frequency_response(G, 10)
print("mag at 10 rad/s:", mag10, "Phase:", ph10)
mag100, ph100, w100 = ctrl.frequency_response(G, 100)
print("mag at 100 rad/s:", mag100, "Phase:", ph100)
[Link](G, dB=True, deg=True,display_margins=True)
[Link]()
[Link]()
gm, pm, wg, wp = [Link](G)
print("Part f - Gain Margin:", gm, "Phase Margin:", pm)
Part a - Closed-loop TF:
<TransferFunction>: sys[79]
Inputs (1): ['u[0]']
Outputs (1): ['y[0]']
25
--------------
s^2 + 5 s + 25
Part b - Poles: [-2.5+4.33012702j -2.5-4.33012702j]
performance metrics {'RiseTime': 0.3349214680718611, 'SettlingTime': 1.6187870956806618,
'SettlingMin': 0.9394891622475168, 'SettlingMax': 1.1630334929041954, 'Overshoot':
16.303349290419543, 'Undershoot': 0.0, 'Peak': 1.1630334929041954, 'PeakTime':
0.7256631808223657, 'SteadyStateValue': 1.0}
mag at 1 rad/s: [4.90290338] Phase: [-1.76819189]
mag at 10 rad/s: [0.2236068] Phase: [-2.67794504]
mag at 100 rad/s: [0.00249688] Phase: [-3.09163426]
Part f - Gain Margin: inf Phase Margin: 51.827292372987756
In [9]:
#question 2
G2 = [Link]([20], [1, 2])
G3 = [Link]([5], [1, 15, 50])
H1 = [Link]([1], [1])
H2 = [Link]([1], [0.1, 1])
H3 = [Link]([0.9], [1])
G1 = [Link]([1], [1])
print("case A:")
Gseries = [Link](G2, G3)
innerloop1 = [Link](Gseries, H1, -1)
innerloop2 = [Link](innerloop1, H2, 1)
outerloop = [Link](innerloop2, H3, -1)
T = [Link](G1, outerloop)
print("- overall tf T(s):", T)
print(" poles:", [Link](T))
print("zeros:", [Link](T))
print(" dc gain:", [Link](T))
t, y = ctrl.step_response(T)
[Link](t, y)
[Link]("Part c: Closed Loop Step Response")
[Link]()
[Link]()
info=ctrl.step_info(T)
print("performance
metrics:",'RiseTime:',info['RiseTime'],',SettlingTime:',info['SettlingTime'],',PeakTime:',info['PeakT
ime'],',Overshoot',info['Overshoot'], ',Peak:',info[
'Peak'],',SteadyStateValue:',info['SteadyStateValue'])
openloop = [Link](G1, innerloop2, H3)
[Link](openloop, dB=True, deg=True,display_margins=True)
[Link]()
[Link]()
gm, pm, wg, wp = [Link](openloop)
print("Gain Margin:", gm, "Phase Margin:", pm)
G1 = [Link]([3.4903047091], [1])
print("case B:")
Gseries = [Link](G2, G3)
innerloop1 = [Link](Gseries, H1, -1)
innerloop2 = [Link](innerloop1, H2, 1)
outerloop = [Link](innerloop2, H3, -1)
T = [Link](G1, outerloop)
print("overall tf T(s):", T)
print(" poles:", [Link](T))
print("zeros:", [Link](T))
print(" dc gain:", [Link](T))
t, y = ctrl.step_response(T)
[Link](t, y)
[Link]("Closed Loop Step Response")
[Link]()
[Link]()
info=ctrl.step_info(T)
print("performance
metrics:",'RiseTime:',info['RiseTime'],',SettlingTime:',info['SettlingTime'],',PeakTime:',info['PeakT
ime'],',Overshoot',info['Overshoot'], ',Peak:',info[
'Peak'],',SteadyStateValue:',info['SteadyStateValue'])
openloop = [Link](G1, innerloop2, H3)
[Link](openloop, dB=True, deg=True,display_margins=True)
[Link]()
[Link]()
gm, pm, wg, wp = [Link](openloop)
print("Gain Margin:", gm, "Phase Margin:", pm)
G1 = [Link]([3, 3], [1, 8])
print("case C:")
Gseries = [Link](G2, G3)
innerloop1 = [Link](Gseries, H1, -1)
innerloop2 = [Link](innerloop1, H2, 1)
outerloop = [Link](innerloop2, H3, -1)
T = [Link](G1, outerloop)
print("overall tf T(s):", T)
print(" poles:", [Link](T))
print("zeros:", [Link](T))
print(" dc gain:", [Link](T))
t, y = ctrl.step_response(T)
[Link](t, y)
[Link]("Closed Loop Step Response")
[Link]()
[Link]()
info=ctrl.step_info(T)
print("performance
metrics:",'RiseTime:',info['RiseTime'],',SettlingTime:',info['SettlingTime'],',PeakTime:',info['PeakT
ime'],',Overshoot',info['Overshoot'], ',Peak:',info[
'Peak'],',SteadyStateValue:',info['SteadyStateValue'])
openloop = [Link](G1, innerloop2, H3)
[Link](openloop, dB=True, deg=True,display_margins=True)
[Link]()
[Link]()
gm, pm, wg, wp = [Link](openloop)
print("Gain Margin:", gm, "Phase Margin:", pm)
case A:
- overall tf T(s): <TransferFunction>: sys[193]
Inputs (1): ['u[0]']
Outputs (1): ['y[0]']
10 s + 100
----------------------------------------
0.1 s^4 + 2.7 s^3 + 25 s^2 + 109 s + 190
poles: [-14.03139728+0.j -4.20314578+3.46609734j
-4.20314578-3.46609734j -4.56231116+0.j ]
zeros: [-10.+0.j]
dc gain: 0.5263157894736842
performance metrics: RiseTime: 0.6142271646225175 ,SettlingTime: 1.0458462532761785
,PeakTime: 1.4276631393928785 ,Overshoot 0.14453911133792352 ,Peak:
0.5270765216386206 ,SteadyStateValue: 0.5263157894736842
Gain Margin: 14.50706493139726 Phase Margin: inf
case B:
overall tf T(s): <TransferFunction>: sys[207]
Inputs (1): ['u[0]']
Outputs (1): ['y[0]']
34.9 s + 349
----------------------------------------
0.1 s^4 + 2.7 s^3 + 25 s^2 + 109 s + 190
poles: [-14.03139728+0.j -4.20314578+3.46609734j
-4.20314578-3.46609734j -4.56231116+0.j ]
zeros: [-10.+0.j]
dc gain: 1.8370024784736843
performance metrics: RiseTime: 0.6142271646225175 ,SettlingTime: 1.0458462532761785
,PeakTime: 1.4276631393928785 ,Overshoot 0.14453911133792122 ,Peak:
1.8396576655313257 ,SteadyStateValue: 1.8370024784736843
Gain Margin: 4.156389238330431 Phase Margin: 64.40923030303077
case C:
overall tf T(s): <TransferFunction>: sys[221]
Inputs (1): ['u[0]']
Outputs (1): ['y[0]']
30 s^2 + 330 s + 300
------------------------------------------------------
0.1 s^5 + 3.5 s^4 + 46.6 s^3 + 309 s^2 + 1062 s + 1520
poles: [-14.03139728+0.j -8. +0.j
-4.20314578+3.46609734j -4.20314578-3.46609734j
-4.56231116+0.j ]
zeros: [-10.+0.j -1.+0.j]
dc gain: 0.19736842105263158
performance metrics: RiseTime: 0.14940660761088287 ,SettlingTime: 1.5106668102878158
,PeakTime: 0.5644249620855576 ,Overshoot 95.22447174128196 ,Peak:
0.38531145738410916 ,SteadyStateValue: 0.19736842105263158
Gain Margin: 14.678901921284181 Phase Margin: inf
In [14]:
#question 3
s = [Link].s
alpha = 86
G = alpha / (s**3 + 14*s**2 + 50*s + 100)
[Link](figsize=(8, 6))
ctrl.root_locus(G)
[Link](f"Root Locus (alpha = {alpha})")
[Link](True)
[Link]()
K1 = 299.71
K2 = 599.42
T1 = [Link](K1 * G, 1)
T2 = [Link](K2 * G, 1)
print(f"Closed-loop poles (K={K1}):\n", [Link](T1))
print(f"\nClosed-loop poles (K={K2}):\n", [Link](T2))
[Link](figsize=(8, 5))
t1, y1 = ctrl.step_response(T1)
t2, y2 = ctrl.step_response(T2)
[Link](t1, y1, label=f"K={K1}")
[Link](t2, y2, label=f"K={K2}")
[Link]("Time (s)")
[Link]("Output")
[Link]("Step Response Comparison")
[Link]()
[Link](True)
[Link]()
[Link]()
ctrl.bode_plot(G, dB=True, margins=True)
[Link](f"Bode Plot (Open Loop, alpha={alpha})")
gm, pm, wgc, wpc = [Link](G)
print("-" * 30)
print("MARGIN ANALYSIS:")
if gm == [Link] or gm == 0:
print("Gain Margin: ∞")
else:
print(f"Gain Margin: {20*np.log10(gm):.2f} dB")
print(f"Phase Margin: {pm:.2f} deg")
print(f"Gain Crossover (wgc): {wgc:.4f} rad/s")
print(f"Phase Crossover (wpc): {wpc:.4f} rad/s")
print("-" * 30)
[Link]()
Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.
Closed-loop poles (K=299.71):
[-34.40550397 +0.j 10.20275198+25.45516646j
10.20275198-25.45516646j]
Closed-loop poles (K=599.42):
[-42.03791621 +0.j 14.0189581 +32.12669273j
14.0189581 -32.12669273j]
C:\Users\1807a\AppData\Local\Temp\ipykernel_28068\[Link]: FutureWarning:
keyword 'margins' is deprecated; use 'display_margins'
ctrl.bode_plot(G, dB=True, margins=True)
------------------------------
MARGIN ANALYSIS:
Gain Margin: 16.87 dB
Phase Margin: inf deg
Gain Crossover (wgc): 7.0711 rad/s
Phase Crossover (wpc): nan rad/s
------------------------------
In [18]:
#question 4
tau = 23
I2 = 9631
K_val = 408
wheel_fwd = [Link]([1], [tau])
wheel_fb = [Link]([1], [1, 0])
wheel_tf = [Link](wheel_fwd, wheel_fb, -1)
pitch_dyn = [Link]([1], [I2, 0, 0])
comp_nogain = [Link]([1, 0.01], [1, 0])
G_ol = [Link](comp_nogain, wheel_tf, pitch_dyn)
[Link]()
ctrl.root_locus(G_ol)
[Link]()
[Link]()
G_actual = [Link]([Link]([K_val], [1]), G_ol)
T_actual = [Link](G_actual, 1, -1)
p = [Link](T_actual)
print("Poles:", p)
p_complex = p[[Link]([Link](p)) > 1e-3]
p1 = p_complex[0]
wn = [Link](p1)
zeta = -[Link](p1) / wn
approx_num = [wn**2]
approx_den = [1, 2 * zeta * wn, wn**2]
T_approx = [Link](approx_num, approx_den)
t_act, y_act = ctrl.step_response(T_actual)
t_app, y_app = ctrl.step_response(T_approx)
[Link]()
[Link](t_act, y_act, label="Actual System")
[Link](t_app, y_app, label="2nd Order Approx", linestyle="--")
[Link]("Step Response Comparison")
[Link]('Time (s)')
[Link]('Output')
[Link]()
[Link]()
[Link]()
Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
Poles: [-0.01539275+0.03484521j -0.01539275-0.03484521j -0.01269276+0.j
0. +0.j ]
In [ ]: