<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Setup Jupyter Notebook
Step1:  cf. Discrete control #1
Step2:  This is hard to see, let's compare to just a ramp
Step3:  Steady State Error is 0.2 (20%), greater than say requirement of Steady State Error being < 2% for ramp input.
Since it's first-order, open loop, it has infinite phase margin, no amount of phase lag will make system unstable.
Create the controller transfer function.
Step4:  Create the feedback system and find the closed loop transfer function.
Step5:  Check the ramp steady state error like we did before
Step6:  There's a jump in error in the beginning, and then goes to 2 percent error.
To find phase margin we look at Bode plot of open loop system.
Step7:  We can convert to a discrete transfer function in the z domain using c2d() 
Default conversion method is Zero Order Hold (ZOH)
Step8:  cf. Discrete control #2
Step9:  Discretization is done with the function c2d
Step10:  cf. Discrete control #3
Step11:  Discrete control #4
Step12:  cf. Discrete control #5
Step13:  cf. Discrete control #6
Step14:  Notch at 0.1 radians per second.
Step15:  Blue digital filter lies on top of analog filter, good.
Change W0, critical frequency.
<ASSISTANT_TASK:>
Python Code:
from pathlib import Path
import sys
notebook_directory_parent = Path.cwd().resolve().parent.parent
if str(notebook_directory_parent) not in sys.path:
    sys.path.append(str(notebook_directory_parent))
print(notebook_directory_parent)
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import control
import control.matlab as controlmatlab
Explanation: Setup Jupyter Notebook
End of explanation
G = control.tf([1,], [0.2, 1])
print(G)
t, y = control.step_response(G)
plt.plot(t, y)
plt.xlabel('Time')
plt.title('Step Response')
unit_step = control.tf([1,], [1, 0])
print(unit_step)
print(G * unit_step)
t, y = control.step_response(G * unit_step)
plt.plot(t, y)
plt.xlabel('Time')
plt.title('Ramp Response')
Explanation: cf. Discrete control #1: Introduction and overview, Aug. 11, 2017, Brian Douglas
End of explanation
t_diff, y_diff = control.step_response(unit_step - G * unit_step)
plt.plot(t_diff, y_diff)
plt.xlabel('Time')
plt.title('Ramp Response Difference')
Explanation: This is hard to see, let's compare to just a ramp
End of explanation
C = control.tf([500, 50], [100, 1, 0])
print(C)
Explanation: Steady State Error is 0.2 (20%), greater than say requirement of Steady State Error being < 2% for ramp input.
Since it's first-order, open loop, it has infinite phase margin, no amount of phase lag will make system unstable.
Create the controller transfer function.
End of explanation
CL = control.feedback(C*G, 1)
print(CL)
Explanation: Create the feedback system and find the closed loop transfer function.
End of explanation
t, y = control.step_response(unit_step - CL * unit_step)
plt.plot(t, y)
plt.xlabel('Time')
plt.title('ramp steady state error')
Explanation: Check the ramp steady state error like we did before
End of explanation
mag, phase, omega = control.matlab.bode(C*G)
# wpc = crossover frequency, wgc = gain crossover frequency (where gain crosses 1)
gain_margin, phase_margin, wpc, wgc = control.margin(C*G)
print(gain_margin)
print(phase_margin) # 50.5 phase margin, beating a 48 degree phase margin
print(wpc)
print(wgc)
Explanation: There's a jump in error in the beginning, and then goes to 2 percent error.
To find phase margin we look at Bode plot of open loop system.
End of explanation
# c2d returns sysd - Discrete time system, with sampling rate Ts
Cz = control.matlab.c2d(C, 0.1) # default 'zoh'
print(Cz)
mag, phase, omega = control.matlab.bode(C, Cz)
gain_margin, phase_margin, wpc, wgc = control.margin(Cz)
print(gain_margin)
print(phase_margin) 
print(wpc)
print(wgc)
Explanation: We can convert to a discrete transfer function in the z domain using c2d() 
Default conversion method is Zero Order Hold (ZOH)
End of explanation
G_s = control.tf([1,], [1, 3, 2])
print(G_s)
Z_g_k = control.matlab.c2d(G_s, 1., 'zoh')
print(Z_g_k)
np.exp(-0.1)
G = control.tf([1], [1, 1])
print(G)
Ts = 0.1 # seconds
Explanation: cf. Discrete control #2: Discretize! Going from continuous to discrete domain Brian Douglas
End of explanation
help(control.matlab.c2d)
Gz = control.matlab.c2d(G, Ts, 'zoh')
print(Gz)
Gf = control.matlab.c2d(G, Ts, 'foh')
print(Gf)
Gi = control.matlab.c2d(G, Ts, 'impulse')
print(Gi)
Gt = control.matlab.c2d(G, Ts, 'tustin')
print(Gt)
Gm = control.matlab.c2d(G, Ts, 'matched')
print(Gm)
Gb = control.matlab.c2d(G, Ts, 'bilinear')
print(Gb)
mag, phase, omega = control.matlab.bode(G, Gz, Gf, Gi, Gt, Gm)
plt.legend()
mag, phase, omega = control.bode_plot([G, Gz, Gf, Gi, Gt, Gm])
plt.gca().legend(('Continuous', 'ZOH', 'FOH', 'Impulse', 'Tustin', 'Matched'))
y, T = control.matlab.step(G)
yz, Tz = control.matlab.step(Gz)
yf, Tf = control.matlab.step(Gf)
yi, Ti = control.matlab.step(Gi)
yt, Tt = control.matlab.step(Gt)
ym, Tm = control.matlab.step(Gm)
plt.plot(T, y, 'r')
plt.plot(Tz, yz, 'b')
plt.plot(Tf, yf, 'g')
plt.plot(Ti, yi, 'y')
plt.plot(Tt, yt, 'p')
plt.plot(Tm, ym, 'o')
plt.show()
print(Gi)
T, y = control.impulse_response(G)
Ti, yi = control.impulse_response(Gi)
Tz, yz = control.impulse_response(Gz)
plt.plot(T, y, 'r')
plt.plot(Tz, yz, 'b')
plt.plot(Ti, yi, 'g')
plt.xlim([0.0, 0.7])
plt.show()
Explanation: Discretization is done with the function c2d
End of explanation
C = control.tf([500, 50], [100, 1, 0])
print(C)
Ct = control.matlab.c2d(C, 0.2, 'tustin')
print(Ct)
Explanation: cf. Discrete control #3: Designing for the zero-order hold Brian Douglas
End of explanation
# No zeroes, poles at -1, -2, gain of 3
b, a = control.matlab.zpk2tf([], [-1, -2], 3)
print(b)
print(a) # Confirms this is same transfer function as written before.
zero_frequency_DC_gain = control.dcgain(control.tf([3], [1,3,2]))
print(zero_frequency_DC_gain)
Gz = control.TransferFunction([1], [1, -0.5032, 0.04978], 1)
print(Gz)
G = control.tf([3,], [1, 3, 2])
print(G)
mag, phase, omega = control.matlab.bode(G, Gz)
t, y = control.step_response(G)
plt.plot(t, y)
plt.xlabel('Time')
plt.title('Step Response')
t, y = control.step_response(Gz)
plt.plot(t, y)
plt.xlabel('Time')
plt.title('Step Response')
t, y = control.step_response(G)
tz, yz = control.step_response(Gz)
plt.plot(t, y)
plt.step(tz, yz)
plt.xlabel('Time')
plt.title('Step Response')
Gz_proper = control.tf([1, 1], [1, -0.5032, -.04978], 1)
mag, phase, omega = control.matlab.bode(G, Gz_proper)
mag, phase, omega = control.matlab.bode(G, Gz, Gz_proper)
t, y = control.step_response(G)
tz, yz = control.step_response(Gz)
tzp, yzp = control.step_response(Gz_proper)
plt.plot(t, y)
plt.step(tz, yz)
plt.step(tzp, yzp)
plt.xlabel('Time')
plt.title('Step Response')
import sympy
z = sympy.Symbol('z')
print( ((z - 0.3679) * (z -0.1353)).as_poly() )
print(G)
print(G.zero())
print(G.pole())
Gz_matched = control.matlab.c2d(G, 1.0, 'matched')
print(Gz_matched)
mag, phase, omega = control.matlab.bode(G, Gz_matched)
t, y = control.step_response(G)
tm, ym = control.step_response(Gz_matched)
plt.plot(t, y)
plt.step(tm, ym)
plt.xlabel('Time')
plt.title('Step Response')
Gz_matched_hand = control.tf([0.4099, 0.4099], [1, -0.5032, 0.04977], 1)
print(Gz_matched_hand)
mag, phase, omega = control.matlab.bode(G, Gz_matched_hand)
t, y = control.step_response(G)
tm, ym = control.step_response(Gz_matched_hand)
plt.plot(t, y)
plt.step(tm, ym)
plt.xlabel('Time')
plt.title('Step Response')
Explanation: Discrete control #4: Discretize with the matched method Brian Douglas
End of explanation
G = control.tf([2,], [1, 1])
print(G)
Gz = control.matlab.c2d(G, 1, 'tustin')
print(Gz)
Gz = control.matlab.c2d(G, 0.1, 'tustin')
print(Gz)
mag, phase, omega = control.matlab.bode(G, Gz)
t, y = control.step_response(G)
tm, ym = control.step_response(Gz)
plt.plot(t, y)
plt.step(tz, yz)
plt.xlabel('Time')
plt.title('Step Response')
Explanation: cf. Discrete control #5: The bilinear transform
End of explanation
G = control.tf([1, 0, 0.01], [1, 0.1, 0.01])
print(G)
# Set sample time to 2 seconds
T = 2
# Set critical frequency to 0.1 rad/s
W0 = 0.1
# Set quality factor to 1
Q = 1
# Create our s-domain notch filter
G1 = control.tf([1, 0, W0*W0], [1, W0/Q, W0*W0])
print(G1)
mag, phase, omega = control.matlab.bode(G1)
Explanation: cf. Discrete control #6: z-plane warping and the bilinear transform
End of explanation
# Let's use the bilinear transform (Tustin) to convert to z-domain
G1t = control.matlab.c2d(G1, T, 'tustin')
print(G1t)
mag, phase, omega = control.matlab.bode(G1, G1t)
Explanation: Notch at 0.1 radians per second.
End of explanation
W0 = 0.7
G1 = control.tf([1, 0, W0*W0], [1, W0/Q, W0*W0])
print(G1)
G1t = control.matlab.c2d(G1, T, 'tustin')
print(G1t)
mag, phase, omega = control.matlab.bode(G1, G1t)
Wa = 0.1 # rad / sec
Wd = 1/(2j) * np.log((1 + Wa * 1j) / (1 - Wa * 1j))
print(Wd)
Wa = 0.7 # rad / sec
Wd = 1/(2j) * np.log((1 + Wa * 1j) / (1 - Wa * 1j))
print(Wd)
Wa = 2.5 # rad / sec
Wd = 1/(2j) * np.log((1 + Wa * 1j) / (1 - Wa * 1j))
print(Wd)
W0 = 0.7
T = 2
Q = 1
G7 = control.tf([1, 0, W0*W0], [1, W0/Q, W0*W0])
print(G7)
# G7 is our s-domain notch filter at 0.7 rad/s
# Now we'll convert it to z-domain with bilinear transform and prewarping
#G7t = control.matlab.c2d(G7, T, ['Method', 'tustin', 'PrewarpFrequency',W0])
G7t = control.sample_system(G7, T, 'tustin', prewarp_frequency=W0)
print(G7t)
mag, phase, omega = control.matlab.bode(G7, G7t)
Explanation: Blue digital filter lies on top of analog filter, good.
Change W0, critical frequency.
End of explanation
<END_TASK>