Scribd
Upload
52 views15 pages

Midterm 2 Codes

The document discusses installing Python packages for control systems analysis and designing a lag-lead compensator for a given system. Specifically, it: 1) Installs the 'control' and 'matlab' Python packages. 2) Calculates damping ratio and natural frequency values to meet design specifications for overshoot, settling time, and steady-state error. 3) Implements a lag-lead compensator with placed poles and zeros, and adjusts the system gain. 4) Plots the root locus, Bode plot, and step response of the compensated system to analyze the design.

Uploaded by

sameertarda
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
52 views15 pages

Midterm 2 Codes

The document discusses installing Python packages for control systems analysis and designing a lag-lead compensator for a given system. Specifically, it: 1) Installs the 'control' and 'matlab' Python packages. 2) Calculates damping ratio and natural frequency values to meet design specifications for overshoot, settling time, and steady-state error. 3) Implements a lag-lead compensator with placed poles and zeros, and adjusts the system gain. 4) Plots the root locus, Bode plot, and step response of the compensated system to analyze the design.

Uploaded by

sameertarda
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd

midterm-2-codes

December 2, 2023

0.1 Control Systems Design (EENG5310-001 & 605)


0.1.1 Take-Home Midterm Examination–2 (Fall 2023)

[3]: !pip install control

Collecting control
Downloading [Link] (455 kB)
���������������������������������������� 455.1/455.1
kB 6.8 MB/s eta [Link]
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-
packages (from control) (1.23.5)
Requirement already satisfied: scipy>=1.3 in /usr/local/lib/python3.10/dist-
packages (from control) (1.11.3)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-
packages (from control) (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in
/usr/local/lib/python3.10/dist-packages (from matplotlib->control) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-
packages (from matplotlib->control) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in
/usr/local/lib/python3.10/dist-packages (from matplotlib->control) (4.44.3)
Requirement already satisfied: kiwisolver>=1.0.1 in
/usr/local/lib/python3.10/dist-packages (from matplotlib->control) (1.4.5)
Requirement already satisfied: packaging>=20.0 in
/usr/local/lib/python3.10/dist-packages (from matplotlib->control) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-
packages (from matplotlib->control) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in
/usr/local/lib/python3.10/dist-packages (from matplotlib->control) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in
/usr/local/lib/python3.10/dist-packages (from matplotlib->control) (2.8.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-
packages (from python-dateutil>=2.7->matplotlib->control) (1.16.0)
Installing collected packages: control
Successfully installed control-0.9.4

[4]: !pip install matlab

1
Collecting matlab
Downloading [Link] (538 bytes)
Preparing metadata ([Link]) … done
Building wheels for collected packages: matlab
Building wheel for matlab ([Link]) … done
Created wheel for matlab: filename=[Link] size=1156
sha256=e25a1bdd5ec7b87af81ca2a8cef0cd4ba7f7bcd98e4935bb9616e48e0b799a30
Stored in directory: /root/.cache/pip/wheels/d1/d3/76/5314058ee22e7957a18eec02
91788462f1df178bb08223bdd2
Successfully built matlab
Installing collected packages: matlab
Successfully installed matlab-0.1

1 QUESTION 1
[54]: from [Link] import *
import [Link] as plt
import control

# Define the transfer function G(s) = K / [s(s+6)(s+10)]


# Since we are plotting the root locus, we can just set K = 1 for now
s = [Link]('s')
G_s = 1 / (s * (s + 6) * (s + 10))

# Plot the root locus


#rlocus = control.root_locus(G_s)

# Plotting the root locus


[Link](figsize=(10, 6))
rlocus(G_s)
[Link]('Root Locus of Given System G(s) = 10/((s+6)*(s+1))\n')
[Link]('Real Axis')
[Link]('Imaginary Axis')
[Link](True)

# Show the plot


[Link]()

# Show the plot


[Link]()

2
[55]: from [Link] import fsolve
import numpy as np

# Define the equation to solve for the damping ratio (zeta) given the overshoot␣
↪(OS)

def overshoot_to_damping_ratio(zeta, OS=0.20):


return [Link](-zeta * [Link] / [Link](1 - zeta**2)) - OS

# Use fsolve to find the damping ratio that corresponds to 20% overshoot
zeta_initial_guess = 0.5
damping_ratio = fsolve(overshoot_to_damping_ratio, zeta_initial_guess)[0]
print(f"The damping ratio for 20% overshoot is: {damping_ratio}")

The damping ratio for 20% overshoot is: 0.4559498107691263

[63]: import numpy as np


from scipy import signal
from [Link] import fsolve
from [Link] import place_poles
import [Link] as plt
from [Link] import *
from control import root_locus
from control import bode

3
from control import step_response

# Define the transfer function G(s) = K / [s(s+6)(s+10)]


# Since we are plotting the root locus, we can just set K = 1 for now
s = [Link]('s')
G_s = 1 / (s * (s + 6) * (s + 10))

# Given system parameters


desired_overshoot = 0.20 # 20% overshoot
desired_settling_time_reduction = 2 # Twofold reduction in settling time
desired_steady_state_error_reduction = 10 # Tenfold reduction in steady-state␣
↪error

# Step 1: Calculate the damping ratio for the desired overshoot


def calculate_damping_ratio(OS):
return [Link]([Link](OS)**2 / ([Link]**2 + [Link](OS)**2))

# Step 2: Calculate the desired natural frequency based on the desired␣


↪reduction in settling time

def calculate_natural_frequency(zeta, Ts_reduction, Ts_original=1):


# Assuming original Ts = 4 / (zeta * wn) for some wn, we get the new wn for␣
↪the reduced Ts.

return 4 / (zeta * Ts_original * Ts_reduction)

# Step 3: Design the lead compensator


# Place the zero at the open-loop pole of -6 (as per the simplification given)
# The pole needs to be calculated based on the root locus angle criterion,␣
↪which is not trivial without plotting.

# For the purpose of this example, let's assume we've calculated it to be at␣
↪-10 (this would normally be done

# iteratively by checking the root locus plot).


lead_zero = -6
lead_pole = -10

# Step 4: Design the lag compensator


# Place the lag compensator pole at 0.01 (as per the instructions)
# The zero of the lag compensator will be placed at a value that achieves the␣
↪desired steady-state error reduction.

# This value should be calculated based on the specific system response and is␣
↪typically done iteratively.

lag_pole = 0.01
lag_zero = 0.001 # This value is an assumption and should be fine-tuned based␣
↪on the system response

# Step 5: Calculate the steady-state error for a ramp input and adjust for the␣
↪desired reduction

4
def adjust_gain_for_steady_state_error(Kv_reduction, G_s):
# Calculate the necessary Kv
Kv_original = (G_s.horner(0))[0][0].real
desired_Kv = Kv_original * Kv_reduction
# Adjust the gain K to achieve the desired Kv
K_adjusted = min(desired_Kv*60, 1000) # because Kv = K/60 for the given␣
↪G(s)

return K_adjusted

# Implement the compensator design


zeta = calculate_damping_ratio(desired_overshoot)
wn = calculate_natural_frequency(zeta, desired_settling_time_reduction)

# Define the open-loop transfer function G(s) without the gain K


numerator = [1]
denominator = [1, 16, 60, 0]
G_s = tf(numerator, denominator)

# Adjust the gain for the steady-state error


K_adjusted =␣
↪adjust_gain_for_steady_state_error(desired_steady_state_error_reduction, G_s)

# Create the transfer function for the lag-lead compensator


# The compensator transfer function C(s) = (s - zero_lead) / (s - pole_lead) *␣
↪(s - zero_lag) / (s - pole_lag)

C_s = tf([1, -lead_zero], [1, -lead_pole]) * tf([1, -lag_zero], [1, -lag_pole])

# The final open-loop transfer function with compensator and adjusted gain
G_s_compensated = K_adjusted * C_s * G_s

print(G_s, "\n", K_adjusted,"\n", C_s, G_s_compensated)

print("\n\n")
print("======================================================================================"

# Plotting the root locus, Bode plot, or step response of the compensated␣
↪system to fine-tune the compensator:

# Plot the root locus


[Link]()
root_locus(G_s_compensated)
[Link]('Root Locus of the Compensated System')
[Link]()

print("\n\n")
print("======================================================================================"

# Plot the Bode plot

5
[Link]()
mag, phase, omega = bode(G_s_compensated, dB=True, Hz=True, deg=True, Plot=True)
plt.tight_layout()
[Link]()

print("\n\n")
print("======================================================================================"

# Plot the step response


[Link]()
time, response = step_response(G_s_compensated)
[Link](time, response)
[Link]('Step Response of the Compensated System')
[Link]('Time (seconds)')
[Link]('Response')
[Link](True)
[Link]()

1
-------------------
s^3 + 16 s^2 + 60 s

1000

s^2 + 5.999 s - 0.006


---------------------
s^2 + 9.99 s - 0.1

1000 s^2 + 5999 s - 6


---------------------------------------------
s^5 + 25.99 s^4 + 219.7 s^3 + 597.8 s^2 - 6 s

================================================================================
======

6
================================================================================
======
/usr/local/lib/python3.10/dist-packages/control/[Link]: FutureWarning:
'Plot' keyword is deprecated in bode_plot; use 'plot'
[Link]("'Plot' keyword is deprecated in bode_plot; use 'plot'",

7
================================================================================
======

8
2 QUESTION 2
[28]: from [Link] import step, TransferFunction
import numpy as np
import [Link] as plt

# Given system parameters


K = 144 # System gain
alpha_lead = 0.8987412594051678 # Ratio of the zero to the pole of the␣
↪compensator

T_lead = 0.012850915297885821 # Time constant of the compensator


Kc_final = 10 # Gain of the compensator

# Uncompensated system transfer function G(s)


num_uncompensated = [100 * K]
den_uncompensated = [1, 136, 3600, 0]
system_uncompensated = TransferFunction(num_uncompensated, den_uncompensated)

# Time vector for step response (we choose enough time to see the full response)

9
t = [Link](0, 2, 1000)

# Step response for the uncompensated system


t_uncompensated, response_uncompensated = step(system_uncompensated, T=t)

# Compensated system transfer function G(s)Gc(s)


num_compensated = [Kc_final, Kc_final / T_lead]
den_compensated = [1, 1 / (alpha_lead * T_lead)]
compensator = TransferFunction(num_compensated, den_compensated)
system_compensated = TransferFunction([Link](num_uncompensated,␣
↪num_compensated),

[Link](den_uncompensated,␣
↪den_compensated))

# Step response for the compensated system


t_compensated, response_compensated = step(system_compensated, T=t)

# Plot the step responses


[Link](figsize=(12, 6))
[Link](t_uncompensated, response_uncompensated, label='Uncompensated System')
[Link](t_compensated, response_compensated, label='Compensated System')
[Link]('Step Response of Uncompensated vs. Compensated Systems')
[Link]('Time (seconds)')
[Link]('Response')
[Link](True)
[Link]()
[Link]()

# Calculate the percent overshoot and peak time for the uncompensated system
peak_uncompensated = [Link](response_uncompensated)
percent_overshoot_uncompensated = (peak_uncompensated - 1) * 100
peak_time_uncompensated = t_uncompensated[[Link](response_uncompensated)]

# Calculate the percent overshoot and peak time for the compensated system
peak_compensated = [Link](response_compensated)
percent_overshoot_compensated = (peak_compensated - 1) * 100
peak_time_compensated = t_compensated[[Link](response_compensated)]

# Output the results


print(f"Uncompensated System: Percent Overshoot =␣
↪{percent_overshoot_uncompensated}%, Peak Time = {peak_time_uncompensated}s")

print(f"Compensated System: Percent Overshoot =␣


↪{percent_overshoot_compensated}%, Peak Time = {peak_time_compensated}s")

10
Uncompensated System: Percent Overshoot = 684.8888888888632%, Peak Time = 2.0s
Compensated System: Percent Overshoot = 6958.798296093113%, Peak Time = 2.0s

3 QUESTION 3
[27]: import numpy as np
import [Link] as plt

# Constants for the transfer functions


K1_values = [0.5, 2]
K2 = 1

# Define the range of omega for the Nyquist plot


omega = [Link](-1, 3, 1000)
s = 1j * omega # j*omega represents the imaginary axis

# Function to calculate the transfer function


def transfer_function(K1, K2, s):
return K1 * (1 + K2 * s) / (s * (s - 1))

# Plot the Nyquist plot for each set of K1 values


for K1 in K1_values:
Gs = transfer_function(K1, K2, s)
[Link]()
[Link]([Link], [Link], label='Nyquist Plot for K1={}, K2={}'.format(K1,␣
↪K2))

11
[Link]([Link], -[Link], linestyle='--') # Mirror image for symmetry
[Link]([-1], [0], color='red') # Stability reference point (-1,0)
[Link]('Real Axis')
[Link]('Imaginary Axis')
[Link]('Nyquist Plot for K1={}, K2={}'.format(K1, K2))
[Link]()
[Link]()
[Link]('equal') # To maintain the aspect ratio

[Link]()

12
[47]: import numpy as np
import [Link] as plt

# Constants for the transfer functions


K1_values = [0.5, 2]
K2 = 1

# Define the range of omega for the Nyquist plot


omega = [Link](-1, 3, 1000)
s = 1j * omega # j*omega represents the imaginary axis

# Function to calculate the transfer function


def transfer_function(K1, K2, s):
return K1 * (1 + K2 * s) / (s * (s - 1))

# Plot the Nyquist plot for each set of K1 values


for K1 in K1_values:
Gs = transfer_function(K1, K2, s)
[Link]()

13
[Link]([Link], [Link], label='Nyquist Plot for K1={}, K2={}'.format(K1,␣
↪K2))
[Link]([Link], -[Link], linestyle='--') # Mirror image for symmetry
[Link]([-1], [0], color='red', zorder=3) # Stability reference point␣
↪(-1,0)

[Link]('Real Axis')
[Link]('Imaginary Axis')
[Link]('Nyquist Plot for K1={}, K2={}\n'.format(K1, K2))
[Link]()
[Link](True)
[Link]('equal') # To maintain the aspect ratio
[Link](-4, 4) # Set x-axis limits to zoom in around the critical point
[Link](-4, 4) # Set y-axis limits to have a clear view of the critical␣
↪point

[Link]()

14
15

You might also like