- Notifications
You must be signed in to change notification settings - Fork450
Description
Version of the control package: 0.10.1
I would like to plot a Bode plot of two third-order systems.
One system has three poles and no zeros. The other system also has three poles and two zeros.
The phase response is step-like due to the lack of damping.
I would have expected the following phase response and verified it with Matlab:
System 1:
Starts at -90° due to the pole at 0. Constant response up to the resonance frequency, then jumps to -270°. I also get this response with the control package.
System 2:
Starts at -90° due to the pole at 0. Constant response up to the frequency of the double zero, then jumps to +90°. Constant curve up to the double pole position and jump to -90°. However, with the control package, I get jumps to -270° and then -450°. The -270° is equivalent to the +90°, but the “direction of rotation,” i.e., how you get there, is not correct.
The code documenting the transfer functions and the rest of the code can be found below. Is it possible to achieve the expected behavior?
Matlab output:
Python control output:
importnumpyasnpfromscipyimportsignalimportplotly.graph_objectsasgoimportsympyasspfromplotly.subplotsimportmake_subplotsimportmatplotlib.pyplotaspltimportcontrol# ParametersLafe_val=3.1e-3Lg_val=2e-3Cf_val=3.3e-6# Symbolic Definitions=sp.symbols('s')Lafe,Lg,Cf=sp.symbols('Lafe Lg Cf')G_g=1/ (Lafe*Lg*Cf*s**3+ (Lafe+Lg)*s)G_c= (1+Lg*Cf*s**2)/ (Lafe*Lg*Cf*s**3+ (Lafe+Lg)*s)# Display symbolic transfer functionsprint("Symbolic transfer function Gg(s):")sp.pprint(G_g)print("Symbolic transfer function Gc(s):")sp.pprint(G_c)## Transfer Function 1# Calculate poles and zeros# Extract numerator and denominatornum_g_sym,den_g_sym=sp.fraction(G_g)zeros_g=sp.solve(num_g_sym,s)poles_g=sp.solve(den_g_sym,s)# Substitute parameter valuesG_g_num=G_g.subs({Lafe:Lafe_val,Lg:Lg_val,Cf:Cf_val})zeros_g_num= [z.subs({Lafe:Lafe_val,Lg:Lg_val,Cf:Cf_val})forzinzeros_g]poles_g_num= [p.subs({Lafe:Lafe_val,Lg:Lg_val,Cf:Cf_val})forpinpoles_g]# Display zerosprint("\nZeros of the numerator:")ifzeros_g:forz_sym,z_numinzip(zeros_g,zeros_g_num):print(f" s ={z_sym} ≈{sp.N(z_num,4)}")else:print(" No zeros (numerator is constant).")# Display polesprint("\nPoles of the denominator:")ifpoles_g:forp_sym,p_numinzip(poles_g,poles_g_num):print(f" s ={p_sym} ≈{sp.N(p_num,4)}")else:print(" No poles (denominator is constant).")# Extract numerator and denominatornum_g,den_g=sp.fraction(G_g_num)# Calculate coefficientsnum_coeffs_g=sp.Poly(num_g,s).all_coeffs()den_coeffs_g=sp.Poly(den_g,s).all_coeffs()# Convert to floatnum_coeffs_g= [float(c)forcinnum_coeffs_g]den_coeffs_g= [float(c)forcinden_coeffs_g]# Transfer functionsystem_g=control.TransferFunction(num_coeffs_g,den_coeffs_g)## Transfer Function 2# Calculate poles and zeros# Extract numerator and denominatornum_c_sym,den_c_sym=sp.fraction(G_c)zeros_c=sp.solve(num_c_sym,s)poles_c=sp.solve(den_c_sym,s)# Substitute parameter valuesG_c_num=G_c.subs({Lafe:Lafe_val,Lg:Lg_val,Cf:Cf_val})zeros_c_num= [z.subs({Lafe:Lafe_val,Lg:Lg_val,Cf:Cf_val})forzinzeros_c]poles_c_num= [p.subs({Lafe:Lafe_val,Lg:Lg_val,Cf:Cf_val})forpinpoles_c]# Display zerosprint("\nZeros of the numerator:")ifzeros_c:forz_sym,z_numinzip(zeros_c,zeros_c_num):print(f" s ={z_sym} ≈{sp.N(z_num,4)}")else:print(" No zeros (numerator is constant).")# Display polesprint("\nPoles of the denominator:")ifpoles_c:forp_sym,p_numinzip(poles_c,poles_c_num):print(f" s ={p_sym} ≈{sp.N(p_num,4)}")else:print(" No poles (denominator is constant).")# Extract numerator and denominatornum_c,den_c=sp.fraction(G_c_num)# Calculate coefficientsnum_coeffs_c=sp.Poly(num_c,s).all_coeffs()den_coeffs_c=sp.Poly(den_c,s).all_coeffs()# Convert to floatnum_coeffs_c= [float(c)forcinnum_coeffs_c]den_coeffs_c= [float(c)forcinden_coeffs_c]# Transfer functionsystem_c=control.TransferFunction(num_coeffs_c,den_coeffs_c)# Frequency rangefrequencies_rad_s=np.logspace(2,5,num=500)frequencies_pole=np.linspace(0.99*float(sp.Abs(poles_c_num[1])),1.01*float(sp.Abs(poles_c_num[1])),100)frequencies_zero=np.linspace(0.99*float(sp.Abs(zeros_c_num[1])),1.01*float(sp.Abs(zeros_c_num[1])),100)combined_frequencies=np.sort(np.unique(np.concatenate((frequencies_rad_s,frequencies_pole,frequencies_zero, [float(sp.Abs(poles_c_num[1]))], [float(sp.Abs(zeros_c_num[1]))]))))# Plot Bode diagramcontrol.bode_plot([system_g,system_c],combined_frequencies,dB=True,Hz=False,deg=True)plt.show()