From 460c5ddcb1a9844d0bdd5f2e393221773eb94ea5 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Thu, 2 May 2024 00:05:47 +0330 Subject: [PATCH 1/8] Add files via upload I tryed to write python equivalent codes of two of the simulations, DCmotor.m and DCmotor_transfun.m, the second on is succesful how ever in the first one there is an error: "incompatible variables" on the simulation line that i could not solve I would be happy if anyone give me a hand on debuging:) tnx --- Python Equivalence/Chapter 3/DCmotor.py | 30 +++++++++++++++++++ .../Chapter 3/DCmotor_transfun.py | 26 ++++++++++++++++ 2 files changed, 56 insertions(+) create mode 100644 Python Equivalence/Chapter 3/DCmotor.py create mode 100644 Python Equivalence/Chapter 3/DCmotor_transfun.py diff --git a/Python Equivalence/Chapter 3/DCmotor.py b/Python Equivalence/Chapter 3/DCmotor.py new file mode 100644 index 0000000..0948923 --- /dev/null +++ b/Python Equivalence/Chapter 3/DCmotor.py @@ -0,0 +1,30 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import lsim + +# Define system matrices +A = np.array([[0, 1, 0], [0, 0, 4.438], [0, -12, -24]]) +b1 = np.array([[0], [0], [20]]) +b2 = np.array([[0], [-7.396], [0]]) +B = np.hstack((b1, b2)) +C = np.array([[1, 0, 0], [0, 1, 0]]) +D = 0 + +# Time vector +t = np.arange(0, 4, 0.01) + +# Define the input signal +u1 = 3 - 6 * square(2 * np.pi * 4 * t) + +# Perform simulation +t, y, x = lsim((A, B, C, D), U=u1, T=t) + +# Plot the results +plt.plot(t, x[:, 0], 'k', label='theta') +plt.plot(t, x[:, 1], 'k-.', label='omega') +plt.plot(t, x[:, 2], 'k:', label='i') +plt.grid(True) +plt.xlabel('Time (sec)') +plt.ylabel('State variable') +plt.legend() +plt.show() diff --git a/Python Equivalence/Chapter 3/DCmotor_transfun.py b/Python Equivalence/Chapter 3/DCmotor_transfun.py new file mode 100644 index 0000000..1459703 --- /dev/null +++ b/Python Equivalence/Chapter 3/DCmotor_transfun.py @@ -0,0 +1,26 @@ +import numpy as np +import control as ctrl + + +# Define system matrices +A = np.array([[0, 1, 0], [0, 0, 4.438], [0, -12, -24]]) +b1 = np.array([[0], [0], [20]]) +b2 = np.array([[0], [-7.396], [0]]) +B = np.hstack((b1, b2)) +C = np.array([[1, 0, 0]]) +D = np.array([[0, 0]]) + +# Create state-space system +DCM = ctrl.ss(A, B, C, D) + +# Convert to transfer function +DCM_tf = ctrl.ss2tf(DCM) + +# Convert to zero-pole-gain +DCM_zpk = ctrl.ss2zpk(DCM) + +# Print the results +print("Transfer Function:") +print(DCM_tf) +print("\nZero-Pole-Gain:") +print(DCM_zpk) From ce9f5d7153ed5101f262f682105a2c1933578e91 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Sun, 14 Jul 2024 17:54:35 +0330 Subject: [PATCH 2/8] Add files via upload --- Python Equivalent/Chapter 5_6/exp5_3.py | 31 +++++++++++++++++++++++ Python Equivalent/Chapter 5_6/exp5_4.py | 33 +++++++++++++++++++++++++ 2 files changed, 64 insertions(+) create mode 100644 Python Equivalent/Chapter 5_6/exp5_3.py create mode 100644 Python Equivalent/Chapter 5_6/exp5_4.py diff --git a/Python Equivalent/Chapter 5_6/exp5_3.py b/Python Equivalent/Chapter 5_6/exp5_3.py new file mode 100644 index 0000000..e1885b4 --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp5_3.py @@ -0,0 +1,31 @@ +import numpy as np +A = np.array([[-3/2, 1/2], [1/2, -3/2]]) +B = np.array([[1/2], [1/2]]) +C = np.array([[1, -1]]) +def ctrb(A, B): + n = A.shape[0] + m = B.shape[1] + C = np.hstack([np.linalg.matrix_power(A, i) @ B for i in range(n)]) + return C +def obsv(A, C): + n = A.shape[0] + m = C.shape[1] + O = np.vstack([C @ np.linalg.matrix_power(A, i) for i in range(n)]) + return O +def null_space(A, tol=1e-15): + u, s, vh = np.linalg.svd(A) + null_mask = (s <= tol) + null_space = np.compress(null_mask, vh, axis=0) + return null_space.T +Cn = ctrb(A, B) +print("Controllability Matrix:") +print(Cn) +print("Rank of Controllability Matrix:", np.linalg.matrix_rank(Cn)) +print("Null space of Controllability Matrix:") +print(null_space(Cn)) +On = obsv(A, C) +print("\nObservability Matrix:") +print(On) +print("Rank of Observability Matrix:", np.linalg.matrix_rank(On)) +print("Null space of Observability Matrix:") +print(null_space(On)) diff --git a/Python Equivalent/Chapter 5_6/exp5_4.py b/Python Equivalent/Chapter 5_6/exp5_4.py new file mode 100644 index 0000000..b2d165e --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp5_4.py @@ -0,0 +1,33 @@ +import numpy as np +num = [[1, 2], [-1, 1]] +den = [[1, 1], [1, 2], [1, 1], [1, 3]] +def nested_lists_to_arrays(nested_list): + return [np.array(coeff) for coeff in nested_list] +num_np = nested_lists_to_arrays(num) +den_np = nested_lists_to_arrays(den) +def tf(num, den): + num_poly = np.poly1d(num) + den_poly = np.poly1d(den) + return num_poly, den_poly +num_tf, den_tf = tf(num_np[0], den_np[0]) +print("Transfer Function (num/den):") +print(num_tf) +print("---") +print(den_tf) +A = np.array([[-1]]) +B = np.array([[1], [2]]) +C = np.array([[1, 2]]) +D = np.array([[0]]) +def ss2tf(A, B, C, D): + return np.linalg.inv(s*np.eye(A.shape[0]) - A) @ B + D +s = np.array([[0, 1], [-1, 0]]) +my_sys = np.array([[1/(s+1), 2/(s+2)], [-1/(s+1), 1/(s+3)]]) +A_ss = np.array([[0, 1], [-1, -4]]) +B_ss = np.array([[0], [0], [1]]) +C_ss = np.array([[1, 2]]) +D_ss = np.array([[0]]) +print("\nState-Space Matrices:") +print("A_ss:\n", A_ss) +print("B_ss:\n", B_ss) +print("C_ss:\n", C_ss) +print("D_ss:\n", D_ss) From 2933f66110af4f2e1b7e012787344844c052c639 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Tue, 23 Jul 2024 22:22:55 +0330 Subject: [PATCH 3/8] Add files via upload --- Python Equivalent/Chapter 5_6/exp5_6.py | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 Python Equivalent/Chapter 5_6/exp5_6.py diff --git a/Python Equivalent/Chapter 5_6/exp5_6.py b/Python Equivalent/Chapter 5_6/exp5_6.py new file mode 100644 index 0000000..097f233 --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp5_6.py @@ -0,0 +1,9 @@ +import numpy as np +from scipy.linalg import solve_continuous_lyapunov +A = np.array([[-1, -2], [1, -4]]) +Q = np.eye(2) +P = solve_continuous_lyapunov(A.T, -Q) +det_P = np.linalg.det(P) +print("Matrix P (Lyapunov equation solution):") +print(P) +print("\nDeterminant of P:", det_P) \ No newline at end of file From 748decc36cdd35cb651c6fb40e67a1147b054280 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Thu, 1 Aug 2024 12:50:48 +0330 Subject: [PATCH 4/8] Add files via upload --- Python Equivalent/Chapter 5_6/exp4_3.py | 13 ++++++ Python Equivalent/Chapter 5_6/exp4_9.py | 21 ++++++++++ Python Equivalent/Chapter 5_6/exp5_1.py | 54 +++++++++++++++++++++++++ Python Equivalent/Chapter 5_6/exp5_2.py | 22 ++++++++++ 4 files changed, 110 insertions(+) create mode 100644 Python Equivalent/Chapter 5_6/exp4_3.py create mode 100644 Python Equivalent/Chapter 5_6/exp4_9.py create mode 100644 Python Equivalent/Chapter 5_6/exp5_1.py create mode 100644 Python Equivalent/Chapter 5_6/exp5_2.py diff --git a/Python Equivalent/Chapter 5_6/exp4_3.py b/Python Equivalent/Chapter 5_6/exp4_3.py new file mode 100644 index 0000000..718b611 --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp4_3.py @@ -0,0 +1,13 @@ +import numpy as np +from scipy.linalg import null_space +A = np.array([[-3/2, 1/2], [1/2, -3/2]]) +C = np.array([[1, -1]]) +O = np.vstack([C @ np.linalg.matrix_power(A, i) for i in range(A.shape[0])]) +rank_O = np.linalg.matrix_rank(O) +null_O = null_space(O) +print("Observability matrix O:") +print(O) +print("\nRank of O:") +print(rank_O) +print("\nNull space of O:") +print(null_O) diff --git a/Python Equivalent/Chapter 5_6/exp4_9.py b/Python Equivalent/Chapter 5_6/exp4_9.py new file mode 100644 index 0000000..9811ecf --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp4_9.py @@ -0,0 +1,21 @@ +import numpy as np +from scipy.linalg import null_space +A = np.array([[-3/2, 1/2], + [1/2, -3/2]]) +B = np.array([[1/2], [1/2]]) +n = A.shape[0] +C = np.zeros((n, n)) +Cb = B +for i in range(n): + C[:, i:i+1] = Cb + Cb = A @ Cb +rank_C = np.linalg.matrix_rank(C) +null_C = null_space(C) +print("Controllability matrix C:") +print(C) + +print("\nRank of C:") +print(rank_C) + +print("\nNull space of C:") +print(null_C) diff --git a/Python Equivalent/Chapter 5_6/exp5_1.py b/Python Equivalent/Chapter 5_6/exp5_1.py new file mode 100644 index 0000000..74524a8 --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp5_1.py @@ -0,0 +1,54 @@ +import numpy as np +from scipy.signal import StateSpace, tf2ss +A1 = np.array([[0, 1, 0], + [0, 0, 1], + [-5, -11, -6]]) +B1 = np.array([[0], [0], [1]]) +C1 = np.array([[1, 0, 1]]) +D1 = np.array([[0]]) +sys1 = StateSpace(A1, B1, C1, D1) +A_tf1, B_tf1, C_tf1, D_tf1 = tf2ss([1, 0, 1], [1, 6, 11, 5]) +my_sys1 = StateSpace(A_tf1, B_tf1, C_tf1, D_tf1) +A2 = np.array([[0, 0, -5], + [1, 0, -11], + [0, 1, -6]]) +B2 = np.array([[1], [0], [1]]) +C2 = np.array([[0, 0, 1]]) +D2 = np.array([[0]]) +sys2 = StateSpace(A2, B2, C2, D2) +A_tf2, B_tf2, C_tf2, D_tf2 = tf2ss([0, 0, 1], [1, 6, 11, 5]) +my_sys2 = StateSpace(A_tf2, B_tf2, C_tf2, D_tf2) +A3 = np.array([[0, 1, 0], + [0, 0, 1], + [-5, -11, -6]]) +B3 = np.array([[1], [-6], [26]]) +C3 = np.array([[1, 0, 0]]) +D3 = np.array([[0]]) +sys3 = StateSpace(A3, B3, C3, D3) +A_tf3, B_tf3, C_tf3, D_tf3 = tf2ss([1, 0, 0], [1, 6, 11, 5]) +my_sys3 = StateSpace(A_tf3, B_tf3, C_tf3, D_tf3) +A4 = np.array([[0, 0, -5], + [1, 0, -11], + [0, 1, -6]]) +B4 = np.array([[1], [0], [0]]) +C4 = np.array([[1, -6, 26]]) +D4 = np.array([[0]]) +sys4 = StateSpace(A4, B4, C4, D4) +A_tf4, B_tf4, C_tf4, D_tf4 = tf2ss([1, -6, 26], [1, 6, 11, 5]) +my_sys4 = StateSpace(A_tf4, B_tf4, C_tf4, D_tf4) +print("System 1 - State-Space Representation:") +print(sys1) +print("\nMy System 1 - State-Space Representation:") +print(my_sys1) +print("\nSystem 2 - State-Space Representation:") +print(sys2) +print("\nMy System 2 - State-Space Representation:") +print(my_sys2) +print("\nSystem 3 - State-Space Representation:") +print(sys3) +print("\nMy System 3 - State-Space Representation:") +print(my_sys3) +print("\nSystem 4 - State-Space Representation:") +print(sys4) +print("\nMy System 4 - State-Space Representation:") +print(my_sys4) diff --git a/Python Equivalent/Chapter 5_6/exp5_2.py b/Python Equivalent/Chapter 5_6/exp5_2.py new file mode 100644 index 0000000..7e3d6ee --- /dev/null +++ b/Python Equivalent/Chapter 5_6/exp5_2.py @@ -0,0 +1,22 @@ +import numpy as np +from scipy.signal import StateSpace, ss2tf +A1 = np.array([[-1, 1, 0], + [0, -1, 0], + [0, 0, -2]]) +B1 = np.array([[0], [1], [1]]) +C1 = np.array([[4, -8, 9]]) +D1 = np.array([[0]]) +sys1 = StateSpace(A1, B1, C1, D1) +num1, den1 = ss2tf(A1, B1, C1, D1) +print("System 1 - Transfer Function:") +print(num1, "/", den1) +a2 = np.array([[-1, 0, 0], + [1, -1, 0], + [0, 0, -2]]) +b2 = np.array([[4], [-8], [9]]) +c2 = np.array([[0, 1, 1]]) +d2 = np.array([[0]]) +sys2 = StateSpace(a2, b2, c2, d2) +num2, den2 = ss2tf(a2, b2, c2, d2) +print("\nSystem 2 - Transfer Function:") +print(num2, "/", den2) From c07a653a5aff0d7398922894d730e9e1b80c0395 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Thu, 8 Aug 2024 23:12:15 +0330 Subject: [PATCH 5/8] Create Chapter 4 opening a new file for python codes that belong to chapter 4 of the bood --- Python Equivalence/Chapter 4 | 1 + 1 file changed, 1 insertion(+) create mode 100644 Python Equivalence/Chapter 4 diff --git a/Python Equivalence/Chapter 4 b/Python Equivalence/Chapter 4 new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/Python Equivalence/Chapter 4 @@ -0,0 +1 @@ + From 944eee1bda56233d108a49eeebe1ee1c4f1eda41 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Thu, 8 Aug 2024 23:23:41 +0330 Subject: [PATCH 6/8] Create 6_2.py --- Python Equivalent/Chapter 7/6_2.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 Python Equivalent/Chapter 7/6_2.py diff --git a/Python Equivalent/Chapter 7/6_2.py b/Python Equivalent/Chapter 7/6_2.py new file mode 100644 index 0000000..72dca99 --- /dev/null +++ b/Python Equivalent/Chapter 7/6_2.py @@ -0,0 +1,27 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Jul 9 11:59:31 2024 + +@author: Lenovo +""" + +import numpy as np +from scipy.signal import place_poles + +# Define the system matrices +A = np.array([[0, 1, 0], + [0, 0, 4.438], + [0, -12, -24]]) + +b = np.array([[0], + [0], + [20]]) + +# Desired closed-loop poles +pd = np.array([-24, -3 - 3j, -3 + 3j]) + +# Compute the state feedback gain using pole placement +k = place_poles(A, b, pd) + +print("State feedback gain (k):") +print(k.gain_matrix) # Print the state feedback gain matrix From f26201f04307d7a2c80f2f02b4ceff0915988aaa Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Thu, 8 Aug 2024 23:24:47 +0330 Subject: [PATCH 7/8] Add files via upload --- Python Equivalent/Chapter 7/6_3.py | 30 +++++++++++++++++++++++ Python Equivalent/Chapter 7/6_4.py | 38 ++++++++++++++++++++++++++++++ Python Equivalent/Chapter 7/6_5.py | 32 +++++++++++++++++++++++++ 3 files changed, 100 insertions(+) create mode 100644 Python Equivalent/Chapter 7/6_3.py create mode 100644 Python Equivalent/Chapter 7/6_4.py create mode 100644 Python Equivalent/Chapter 7/6_5.py diff --git a/Python Equivalent/Chapter 7/6_3.py b/Python Equivalent/Chapter 7/6_3.py new file mode 100644 index 0000000..494101f --- /dev/null +++ b/Python Equivalent/Chapter 7/6_3.py @@ -0,0 +1,30 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Jul 9 12:11:09 2024 + +@author: Lenovo +""" + +import numpy as np +import control as ctrl + +# Define the system matrices +A = np.array([[0, 1, 0, 0], + [0, 0, -9.8, 0], + [0, 0, 0, 1], + [0, 0, 19.6, 0]]) + +b = np.array([[0], + [1], + [0], + [-1]]) + +# Define the weighting matrices +Q = np.diag([4, 0, 8.16, 0]) +R = 1 / 400 + +# Compute the LQR gain +k, _, _ = ctrl.lqr(A, b, Q, R) + +print("State feedback gain (k):") +print(k) diff --git a/Python Equivalent/Chapter 7/6_4.py b/Python Equivalent/Chapter 7/6_4.py new file mode 100644 index 0000000..671b70a --- /dev/null +++ b/Python Equivalent/Chapter 7/6_4.py @@ -0,0 +1,38 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Jul 9 12:11:09 2024 + +@author: Lenovo +""" + +import numpy as np + +# Define the system matrices +A = np.array([[0, 1, 0, 0], + [0, 0, -9.8, 0], + [0, 0, 0, 1], + [0, 0, 19.6, 0]]) + +b = np.array([[0], + [1], + [0], + [-1]]) + +# Controllability matrix +C = np.hstack([b, A @ b, A @ A @ b, A @ A @ A @ b]) + +# Given parameters +a = np.array([0, -19.6, 0, 0]) +alpha = np.array([12.86, 63.065, 149.38, 157.0]) + +# Construct Psi matrix +Psi = np.array([[1, a[0], a[1], a[2]], + [0, 1, a[0], a[1]], + [0, 0, 1, a[0]], + [0, 0, 0, 1]]) + +# Compute the state feedback gain k +k = (alpha - a) @ np.linalg.inv(C @ Psi) + +print("State feedback gain (k):") +print(k) diff --git a/Python Equivalent/Chapter 7/6_5.py b/Python Equivalent/Chapter 7/6_5.py new file mode 100644 index 0000000..7dc220f --- /dev/null +++ b/Python Equivalent/Chapter 7/6_5.py @@ -0,0 +1,32 @@ +import numpy as np + +# Define the system matrices +A = np.array([[0, 1, 0, 0], + [0, 0, -9.8, 0], + [0, 0, 0, 1], + [0, 0, 19.6, 0]]) + +b = np.array([[0], + [1], + [0], + [-1]]) + +# Controllability matrix +C = np.hstack([b, A @ b, A @ A @ b, A @ A @ A @ b]) + +# Given parameters +a = np.array([0, -19.6, 0, 0]) +alpha = np.array([12.86, 63.065, 149.38, 157.0]) + +# Construct Psi_1 matrix +Psi_1 = np.array([[1, -a[0], a[0]**2 - a[1], -a[0]**3 + 2*a[0]*a[1] - a[2]], + [0, 1, -a[0], a[0]**2 - a[1]], + [0, 0, 1, -a[0]], + [0, 0, 0, 1]]) + +# Compute the state feedback gain k +k = (alpha - a) @ Psi_1 @ np.linalg.inv(C) + +print("State feedback gain (k):") +print(k) + From dc1d44b530ff76869991ae6efa8daa118ce26855 Mon Sep 17 00:00:00 2001 From: Farbod Raeisi Date: Thu, 8 Aug 2024 23:26:14 +0330 Subject: [PATCH 8/8] Add files via upload --- Python Equivalent/Chapter 7/6_6.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 Python Equivalent/Chapter 7/6_6.py diff --git a/Python Equivalent/Chapter 7/6_6.py b/Python Equivalent/Chapter 7/6_6.py new file mode 100644 index 0000000..e980670 --- /dev/null +++ b/Python Equivalent/Chapter 7/6_6.py @@ -0,0 +1,26 @@ +import numpy as np +import control as ctrl + +# Define the system matrices +A = np.array([[0, 1, 0, 0], + [0, 0, -9.8, 0], + [0, 0, 0, 1], + [0, 0, 19.6, 0]]) + +b = np.array([[0], + [1], + [0], + [-1]]) + +# Calculate eigenvalues of A +e = np.linalg.eigvals(A) +print("Eigenvalues of A:", e) + +# Desired closed-loop poles +pd = np.array([-4.43, -4.43, -2-2j, -2+2j]) + +# Compute the state feedback gain using pole placement (Ackermann method) +k = ctrl.acker(A, b, pd) + +print("State feedback gain (k):") +print(k)