《信号与系统》第八章 系统的状态变量分析

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 电路参数
R = 10  # 电阻(Ω)
L = 0.1  # 电感(H)
C = 1e-3  # 电容(F)

# 状态方程矩阵
A = np.array([
    [0, 1/L],
    [-1/C, -R/L]
])
B = np.array([
    [0],
    [1/L]
])
C = np.array([
    [1, 0]  # 输出为电容电压
])
D = np.array([
    [0]
])

# 定义状态方程的微分形式
def state_eq(t, x):
    # x: 状态变量 [电容电压, 电感电流]
    # u: 输入电压，这里使用外部定义的input_signal函数
    u = input_signal(t)
    dx_dt = A @ x + B.flatten() * u  # 注意这里将B向量展平
    return dx_dt

# 初始状态
x0 = np.array([0, 0])  # 初始电容电压和电感电流为0

# 输入信号：单位阶跃函数
def input_signal(t):
    return 1.0 if t >= 0 else 0.0

# 时间范围
t_span = (0, 1)
t_eval = np.linspace(0, 1, 1000)

# 正确求解状态方程
result = solve_ivp(
    state_eq,
    t_span,
    x0,
    t_eval=t_eval,
    method='RK45'
)

# 提取状态变量
solution = result.y.T  # 转置以获得形状为(n_samples, n_states)的数组

# 计算输出
y = np.array([C @ x for x in solution]).flatten()

# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(t_eval, y, 'b-', linewidth=2)
plt.xlabel('时间 (s)')
plt.ylabel('电容电压 (V)')
plt.title('RLC串联电路阶跃响应')
plt.grid(True)
plt.show()

import sympy as sp

# 定义符号变量
s, t, v_C, i_L, v_in = sp.symbols('s t v_C i_L v_in')
R1, R2, L, C_cap = sp.symbols('R1 R2 L C')  # 将C重命名为C_cap以避免冲突

# 定义表示导数的符号变量
dv_C_dt, di_L_dt = sp.symbols('dv_C_dt di_L_dt')

# 由电路图建立方程，用dv_C_dt和di_L_dt替换导数
eq1 = sp.Eq(C_cap * dv_C_dt, (v_in - v_C)/R1 - i_L)
eq2 = sp.Eq(L * di_L_dt, v_C - i_L * R2)

# 求解状态方程
state_eq = sp.solve([eq1, eq2], [dv_C_dt, di_L_dt])

# 构建状态矩阵A和输入矩阵B
A = sp.Matrix([
    [ -1/(R1*C_cap), -1/C_cap ],
    [ 1/L, -R2/L ]
])
B = sp.Matrix([
    [ 1/(R1*C_cap) ],
    [ 0 ]
])

# 输出方程：假设输出为电感电流
output_matrix = sp.Matrix([[0, 1]])  # 使用行向量表示输出矩阵
D = sp.Matrix([[0]])  # D矩阵也应该是二维的

# 打印状态方程
print("状态方程:")
print(f"dv_C/dt = {state_eq[dv_C_dt]}")
print(f"di_L/dt = {state_eq[di_L_dt]}")

# 打印矩阵形式
print("\n矩阵形式:")
print(f"A = \n{A}")
print(f"B = \n{B}")
print(f"输出矩阵 C = \n{output_matrix}")
print(f"D = \n{D}")

import sympy as sp

# 定义符号变量
s, t = sp.symbols('s t')
u = sp.Function('u')(t)  # u是关于t的函数
y = sp.Function('y')(t)  # y是关于t的函数

# 其中h0是一个常数，选择h0=1（因为输入导数的系数为1）
h0 = 1

# 输入输出微分方程
ode = sp.Eq(sp.diff(y, t, t) + 3*sp.diff(y, t) + 2*y,
            sp.diff(u, t) + 4*u)

# 定义状态变量 - 修正：将h0的定义移到此处之前
x1 = y
x2 = sp.diff(y, t) - h0 * u

# 状态方程
dx1_dt = x2 + h0 * u
dx2_dt = -2*x1 - 3*x2 + (4 - 3*h0) * u

# 构建状态矩阵
A = sp.Matrix([
    [0, 1],
    [-2, -3]
])
B = sp.Matrix([
    [h0],
    [4 - 3*h0]
])
C = sp.Matrix([[1, 0]])  # 输出y = x1
D = sp.Matrix([[h0]])    # D矩阵不为零

# 打印结果
print("状态方程矩阵:")
print(f"A = \n{A}")
print(f"B = \n{B}")
print(f"C = \n{C}")
print(f"D = \n{D}")

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 确保使用交互式后端并指定Qt5Agg（Windows系统推荐）
plt.switch_backend('Qt5Agg')
plt.ion()  # 启用交互模式

# 差分方程系数
a1, a2 = 0.5, 0.06  # y(n) + a1y(n-1) + a2y(n-2) = ...
b0, b1 = 1, 1        # ... = b0u(n) + b1u(n-1)

# 构建正确的状态方程 - 使用一维数组表示行向量
A = np.array([
    [-a1, -a2],
    [1, 0]
])
B = np.array([1, 0])  # 一维数组（行向量）
C = np.array([b1, b0 - a1*b1])  # 一维数组（行向量）
D = np.array([b0 - a1*b1 - a2*b0])  # 一维数组（行向量）

# 系统模拟：计算单位脉冲响应
N = 50  # 计算点数
h = np.zeros(N)
x = np.zeros(2)  # 初始状态

for n in range(N):
    u = 1 if n == 0 else 0  # 单位脉冲输入
    y = C @ x + D * u  # 标量值
    h[n] = y  # 直接赋值
    x = A @ x + B * u  # 更新状态

# 绘制脉冲响应并显示
plt.figure(figsize=(10, 6))
plt.stem(range(N), h, 'b-', linefmt='b-', markerfmt='bo')
plt.xlabel('n')
plt.ylabel('h(n)')
plt.title('离散系统脉冲响应（状态方程法）')
plt.grid(True)
plt.show(block=False)  # 非阻塞显示

# 直接求解差分方程
h_direct = np.zeros(N)
for n in range(N):
    u = 1 if n == 0 else 0
    h_direct[n] = b0*u
    if n >= 1:
        h_direct[n] += b1*u - a1*h_direct[n-1]
    if n >= 2:
        h_direct[n] -= a2*h_direct[n-2]

# 比较两种方法结果
plt.figure(figsize=(10, 6))
plt.stem(range(N), h, 'b-', linefmt='b-', markerfmt='bo', label='状态方程法')
plt.stem(range(N), h_direct, 'r--', linefmt='r--', markerfmt='ro', label='直接法')
plt.xlabel('n')
plt.ylabel('h(n)')
plt.title('脉冲响应对比')
plt.legend()
plt.grid(True)
plt.show(block=False)  # 非阻塞显示

# 计算误差
error = np.sum((h - h_direct)** 2)
print(f"两种方法的均方误差: {error:.6f}")

# 保持程序运行，防止图像窗口关闭
try:
    plt.pause(10)  # 暂停10秒显示图像
except KeyboardInterrupt:
    pass

# 或者使用input()保持程序运行
input("按Enter键退出...")

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import matplotlib as mpl

# 确保使用交互式后端并指定Qt5Agg（Windows系统推荐）
plt.switch_backend('Qt5Agg')
plt.ion()  # 启用交互模式
# 解决中文显示问题和下标问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False         # 解决负号显示问题


# 状态方程矩阵
A = np.array([[-2, 1], [0, -3]])
B = np.array([[1], [1]])
C = np.array([[1, 0]])
D = np.array([[0]])

# 初始状态
x0 = np.array([1, 0])

# 输入信号：指数函数u(t) = e^(-t)u(t)
def u(t):
    return np.exp(-t)

# 使用ODE求解器计算状态响应
def state_equation(t, x):
    return A @ x + B.flatten() * u(t)

# 时间范围
t_span = (0, 5)
t_eval = np.linspace(0, 5, 100)

# 使用ODE求解器求解状态方程
sol = solve_ivp(state_equation, t_span, x0, t_eval=t_eval)
x_states = sol.y.T  # 状态向量，形状为(100, 2)

# 正确计算输出 y(t) = C*x(t) + D*u(t)
u_vals = u(t_eval)
y_output = (C @ x_states.T).flatten() + (D * u_vals).flatten()

# 绘制结果 - 只显示输出
plt.figure(figsize=(10, 6))
plt.plot(t_eval, y_output, 'b-', linewidth=2)
plt.xlabel('时间 (s)', fontsize=12)
plt.ylabel('输出 y(t)', fontsize=12)
plt.title('状态方程求解结果', fontsize=14)
plt.grid(True)
plt.tight_layout()
plt.show()

# 绘制状态变量变化 - 使用普通数字代替下标
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(t_eval, x_states[:, 0], 'r-', label='x1(t)')
plt.plot(t_eval, x_states[:, 1], 'g-', label='x2(t)')
plt.xlabel('时间 (s)', fontsize=12)
plt.ylabel('状态变量', fontsize=12)
plt.legend(fontsize=10)
plt.title('状态变量变化', fontsize=14)
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(t_eval, y_output, 'b-', label='y(t)')
plt.xlabel('时间 (s)', fontsize=12)
plt.ylabel('输出', fontsize=12)
plt.title('系统输出', fontsize=14)
plt.grid(True)
plt.tight_layout()
plt.show()

# 保持程序运行，防止图像窗口关闭
try:
    plt.pause(10)  # 暂停10秒显示图像
except KeyboardInterrupt:
    pass

# 或者使用input()保持程序运行
input("按Enter键退出...")
print("绘图完成！图像已显示")

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
import matplotlib as mpl

# 解决中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 解决负号显示问题
mpl.rcParams['font.family'] = 'sans-serif'    # 设置全局字体

# 状态方程矩阵
A = np.array([
    [-1, 1, 0],
    [0, -2, 1],
    [0, 0, -3]
])
B = np.array([[1], [1], [1]])
C = np.array([[1, 0, 0]])
D = np.array([[0]])

# 计算系统极点（A矩阵特征值）
eigenvalues = linalg.eigvals(A)
real_parts = np.real(eigenvalues)
imag_parts = np.imag(eigenvalues)

# 判断稳定性
is_stable = all(real_parts < 0)

# 绘制极点位置（优化可视化）
plt.figure(figsize=(8, 6))
plt.scatter(real_parts, imag_parts, color='red', s=100, marker='o', zorder=3)

# 绘制坐标轴和网格
plt.axvline(x=0, color='k', linestyle='--', alpha=0.7)  # 虚轴
plt.axhline(y=0, color='k', alpha=0.7)  # 实轴
plt.grid(True, linestyle='--', alpha=0.7)

# 添加极点标签
for i, (re, im) in enumerate(zip(real_parts, imag_parts)):
    if abs(im) < 1e-6:  # 实极点
        plt.text(re+0.1, im+0.1, f'λ{i+1}={re:.1f}', fontsize=10)
    else:  # 复极点
        plt.text(re+0.1, im+0.1, f'λ{i+1}={re:.1f}{"+" if im>=0 else ""}{im:.1f}j', fontsize=10)

# 标记稳定性
stability_text = "系统稳定（所有极点在左半平面）" if is_stable else "系统不稳定（存在右半平面极点）"
plt.text(min(real_parts)-0.5, max(imag_parts)+0.5, stability_text,
         fontsize=12, fontweight='bold',
         bbox=dict(facecolor='white', alpha=0.8))

# 设置坐标轴范围和标签
plt.xlabel('实部', fontsize=12)
plt.ylabel('虚部', fontsize=12)
plt.title('三阶系统极点分布图', fontsize=14, fontweight='bold')
plt.axis('equal')
plt.tight_layout()
plt.show()

# 打印结果（格式化输出）
print("="*50)
print("系统极点分析结果：")
print("-"*50)
for i, val in enumerate(eigenvalues):
    if np.imag(val) == 0:
        print(f"极点 {i+1}: λ={val.real:.2f}（实极点）")
    else:
        print(f"极点 {i+1}: λ={val.real:.2f}{'+' if val.imag>0 else ''}{val.imag:.2f}j（复极点）")
print("-"*50)
print(f"系统稳定性: {'稳定' if is_stable else '不稳定'}")
print("="*50)

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 状态方程矩阵
A = np.array([[0.5, 0.1], [0.2, 0.3]])
B = np.array([1, 0])  # 修改为一维数组
C = np.array([1, 0])  # 修改为一维数组
D = 0  # 标量值

# 初始状态
x = np.array([1.0, 0.0])  # 确保为浮点数


# 输入信号：单位阶跃函数u(n) = u(n)
def u(n):
    return 1.0 if n >= 0 else 0.0  # 返回浮点数


# 计算范围
n_vals = np.arange(0, 20)
y = np.zeros(len(n_vals))

# 存储状态历史
x_history = np.zeros((len(n_vals), 2))

# 计算状态和输出（迭代法）
for i, n_val in enumerate(n_vals):
    # 计算输出 (确保结果为标量)
    y[i] = np.dot(C, x) + D * u(n_val)

    # 保存当前状态
    x_history[i] = x

    # 更新下一状态（如果还没到末尾）
    if i < len(n_vals) - 1:
        # 状态更新: x[n+1] = A @ x[n] + B @ u[n]
        x = np.dot(A, x) + B * u(n_val)

# 绘制结果
plt.figure(figsize=(10, 6))
markerline, stemlines, baseline = plt.stem(n_vals, y, linefmt='b-', markerfmt='bo', basefmt='r-')
plt.setp(markerline, markersize=5)
plt.setp(stemlines, linewidth=1)
plt.xlabel('n')
plt.ylabel('输出 y(n)')
plt.title('离散系统状态方程迭代解法')
plt.grid(True)
plt.show()

# 可选：绘制状态变量变化
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.stem(n_vals, x_history[:, 0], linefmt='r-', markerfmt='ro', label='x1(n)')
plt.stem(n_vals, x_history[:, 1], linefmt='g-', markerfmt='go', label='x2(n)')
plt.xlabel('n')
plt.ylabel('状态变量')
plt.legend()
plt.title('状态变量变化')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.stem(n_vals, y, linefmt='b-', markerfmt='bo', label='y(n)')
plt.xlabel('n')
plt.ylabel('输出')
plt.title('系统输出')
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
import matplotlib as mpl

# 解决中文显示问题
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False    # 解决负号显示问题
mpl.rcParams['font.family'] = 'sans-serif'    # 设置全局字体

# 状态方程矩阵
A = np.array([
    [0.5, 0.3],
    [-0.2, 0.4]
])
B = np.array([[1], [1]])
C = np.array([[1, 0]])
D = np.array([[0]])

# 计算系统函数矩阵H(z)的极点
# 极点为A矩阵的特征值
eigenvalues = linalg.eigvals(A)

# 判断稳定性
is_stable = all(np.abs(eigenvalues) < 1)

# 绘制极点位置（单位圆内）
theta = np.linspace(0, 2*np.pi, 100)
plt.figure(figsize=(8, 8))

# 绘制单位圆
plt.plot(np.cos(theta), np.sin(theta), 'k--', linewidth=1.5, label='单位圆')

# 绘制极点
plt.plot(np.real(eigenvalues), np.imag(eigenvalues), 'ro', markersize=10, label='极点')

# 添加坐标轴
plt.axhline(y=0, color='k', linestyle='-', alpha=0.5)  # 实轴
plt.axvline(x=0, color='k', linestyle='-', alpha=0.5)  # 虚轴

# 设置网格和标签
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('实部')
plt.ylabel('虚部')
plt.title('离散系统极点分布')
plt.axis('equal')

# 设置坐标轴范围
plt.xlim(-1.2, 1.2)
plt.ylim(-1.2, 1.2)

# 标记稳定性
if is_stable:
    plt.text(-1.15, 1.15,
             '系统稳定（所有极点在单位圆内）',
             fontsize=12, color='g')
else:
    plt.text(-1.15, 1.15,
             '系统不稳定（存在单位圆外极点）',
             fontsize=12, color='r')

# 添加图例
plt.legend(loc='best')

plt.tight_layout()
plt.show()

# 打印结果
print("="*50)
print("离散系统稳定性分析结果：")
print("-"*50)
print(f"A矩阵特征值: {eigenvalues}")
print(f"特征值模长: {np.abs(eigenvalues)}")
print(f"系统稳定性: {'稳定' if is_stable else '不稳定'}")
print("="*50)

import numpy as np
import sys
from scipy import linalg

# 状态方程矩阵
A = np.array([
    [0, 1, 0],
    [0, 0, 1],
    [-6, -11, -6]
])
B = np.array([[0], [0], [1]])
C = np.array([[1, 0, 0]])

# 计算可控性矩阵M
n = A.shape[0]  # 状态变量维数
M = np.hstack([np.linalg.matrix_power(A, i) @ B for i in range(n)])

# 计算可观测性矩阵N
N = np.vstack([C @ np.linalg.matrix_power(A, i) for i in range(n)])

# 计算秩
rank_M = np.linalg.matrix_rank(M)
rank_N = np.linalg.matrix_rank(N)

# 判断可控性和可观测性
is_controllable = rank_M == n
is_observable = rank_N == n

# 打印结果
print("=" * 60)
print("系统可控性与可观测性分析")
print("=" * 60)
print("可控性矩阵M:")
print(M)
print(f"可控性矩阵秩: {rank_M} (系统维度: {n})")
print(f"系统可控性: {'可控' if is_controllable else '不可控'}")

print("\n可观测性矩阵N:")
print(N)
print(f"可观测性矩阵秩: {rank_N} (系统维度: {n})")
print(f"系统可观测性: {'可观测' if is_observable else '不可观测'}")
print("=" * 60)

# 线性变换示例：对角化A矩阵
# 计算特征值和特征向量
eigenvalues, V = linalg.eig(A)

# 检查特征向量矩阵是否可逆
if np.linalg.cond(V) < 1 / sys.float_info.epsilon:
    # 变换矩阵P（特征向量矩阵）
    P = V

    # 变换后矩阵
    A_bar = np.linalg.inv(P) @ A @ P
    B_bar = np.linalg.inv(P) @ B
    C_bar = C @ P

    print("\n线性变换后矩阵:")
    print("A_bar =")
    print(np.real_if_close(A_bar))  # 只显示实数部分（如果接近实数）
    print("\nB_bar =")
    print(np.real_if_close(B_bar))
    print("\nC_bar =")
    print(np.real_if_close(C_bar))

    # 检查变换后系统特性是否保持
    print("\n变换后系统特性验证:")
    M_bar = np.hstack([np.linalg.matrix_power(A_bar, i) @ B_bar for i in range(n)])
    N_bar = np.vstack([C_bar @ np.linalg.matrix_power(A_bar, i) for i in range(n)])

    print(f"变换后可控性矩阵秩: {np.linalg.matrix_rank(M_bar)}")
    print(f"变换后可观测性矩阵秩: {np.linalg.matrix_rank(N_bar)}")
    print("=" * 60)
else:
    print("\n警告：特征向量矩阵不可逆，无法进行对角化变换")
    print("=" * 60)

# 系统实现分析
print("\n系统实现分析:")
if is_controllable and is_observable:
    print("系统是既可控又可观测的最小实现")
elif is_controllable:
    print("系统可控但不可观测")
elif is_observable:
    print("系统可观测但不可控")
else:
    print("系统既不可控也不可观测")
print("=" * 60)