模糊PID控制加热器系统仿真

jmfans 发布于 2024-10-12 527 次阅读

AI 摘要

加热器系统仿真采用模糊PID控制以优化温度调节，减少误差。

模糊PID控制加热器系统仿真

目录

引言

本文提供了由自适应模糊PID控制器和Smith预估器控制的加热器系统仿真模型的完整算法仿真。该系统旨在通过根据设定值变化调整加热器的功率输出，最小化超调量、稳态误差和稳定时间来调节加热器的温度。

系统概述

该仿真包含几个相互连接的组件：

  1. 加热器系统: 模拟加热器的热力学特性，包括热输入、环境热损失和功率限制。
  2. PID控制器: 传统的比例-积分-微分控制器，根据误差反馈计算控制信号。
  3. 模糊PID控制器: 通过使用模糊逻辑根据当前误差和误差率自适应调整PID参数（Kp、Ki、Kd）来增强PID控制器。
  4. Smith预估器: 补偿系统延迟，并在存在死区时间的情况下提高控制器的性能。
  5. 仿真框架: 协调这些组件之间的交互，应用设定值更改并记录系统性能指标。
  6. 用户界面: 可视化系统响应、控制信号、PID参数和性能指标，并允许通过交互式滑块进行实时调整。

数学模型

隶属度函数

隶属度函数（MFs）是模糊逻辑的基础，定义了每个输入在每个模糊集中的隶属程度。该仿真使用三种类型的隶属度函数：S形、Z形和三角形。

S形函数

定义:

S形函数用于模拟从低隶属度值逐渐增加到高隶属度值的模糊集。

对于给定的输入 $x$，参数 $a$ 和 $b$ 定义了过渡：

$$
\text{s_shaped}(x, a, b) =
\begin{cases}
0 & \text{if } x \leq a \
2 \left( \frac{x - a}{b - a} \right)^2 & \text{if } a < x \leq \frac{a + b}{2} \
1 - 2 \left( \frac{b - x}{b - a} \right)^2 & \text{if } \frac{a + b}{2} < x < b \
1 & \text{if } x \geq b
\end{cases}
$$

Z形函数

定义:

Z形函数模拟从高隶属度值平滑递减到低隶属度值的模糊集。

对于给定的输入 $x$，参数 $a$ 和 $b$ 定义了过渡：

$$
\text{z_shaped}(x, a, b) =
\begin{cases}
1 & \text{if } x \leq a \
1 - 2 \left( \frac{x - a}{b - a} \right)^2 & \text{if } a < x \leq \frac{a + b}{2} \
2 \left( \frac{b - x}{b - a} \right)^2 & \text{if } \frac{a + b}{2} < x < b \
0 & \text{if } x \geq b
\end{cases}
$$

三角形函数

定义:

三角形函数创建一个简单的、对称的模糊集，在中心点达到峰值。

对于给定的输入 $x$，参数 $a$, $b$, 和 $c$ 定义了三角形：

$$
\text{tri}(x, a, b, c) =
\begin{cases}
0 & \text{if } x \leq a \text{ or } x \geq c \
\frac{x - a}{b - a} & \text{if } a < x \leq b \
\frac{c - x}{c - b} & \text{if } b < x < c \
0 & \text{otherwise}
\end{cases}
$$

模糊PID控制器

模糊PID控制器使用模糊逻辑根据当前误差 ($E$) 和误差导数 ($dE$) 自适应地调整PID参数。

模糊化

过程:

  1. 输入变量:

    • 误差 ($E$): 设定值与过程变量之间的偏差。
    • 误差导数 ($dE$): 误差的变化率。

  2. 模糊化:

    • 每个输入变量都被分配到七个模糊集的隶属度：NB（负大）、NM（负中）、NS（负小）、ZE（零）、PS（正小）、PM（正中）、PB（正大）。
    • 使用上面定义的隶属度函数计算隶属度。

数学表示:

$$
\muE(\text{Label}) = \text{MembershipFunction}(E, \text{parameters})
$$
$$
\mu{dE}(\text{Label}) = \text{MembershipFunction}(dE, \text{parameters})
$$

规则库

结构:

一个 7x7 规则矩阵定义了 $E$ 和 $dE$ 的组合如何映射到 PID 参数的调整 ($\Delta K_p$, $\Delta K_i$, $\Delta K_d$)。

规则示例:

如果 $E$ 为 NB 且 $dE$ 为 PB，则：

$$
\Delta K_p = \text{output_map_Kp}['PB'] = 0.5
$$
$$
\Delta K_i = \text{output_map_Ki}['ZE'] = 0.0
$$
$$
\Delta K_d = \text{output_map_Kd}['NS'] = -0.05
$$

规则应用:

对于 $E$ 和 $dE$ 的每个组合：

  1. 确定隶属度 $\muE[\text{Label}]$ 和 $\mu{dE}[\text{Label}]$。
  2. 应用最小运算符计算触发强度 $\alpha = \min(\mu_E[\text{Label}E], \mu{dE}[\text{Label}_{dE}])$。
  3. 将对 $\Delta K_p$, $\Delta K_i$, 和 $\Delta K_d$ 的贡献加权 $\alpha$ 并聚合。

解模糊

方法: 加权平均

对于每个 PID 参数调整：

$$
\Delta K_p = \frac{\sum (\Delta K_p^{(i)} \cdot \alpha_i)}{\sum \alpha_i}
$$
$$
\Delta K_i = \frac{\sum (\Delta K_i^{(i)} \cdot \alpha_i)}{\sum \alpha_i}
$$
$$
\Delta K_d = \frac{\sum (\Delta K_d^{(i)} \cdot \alpha_i)}{\sum \alpha_i}
$$

如果分母为零（没有触发规则），则调整默认为零。

PID控制器

传统的PID控制器根据误差反馈计算控制信号。

控制律:

$$
u(t) = K_p \cdot e(t) + Ki \cdot \int{0}^{t} e(\tau) d\tau + K_d \cdot \frac{de(t)}{dt}
$$

离散实现:

给定采样时间为 $\Delta t$ 的离散时间步长：

$$
\text{积分项} += e(t) \cdot \Delta t
$$
$$
\text{微分项} = \frac{e(t) - e(t-1)}{\Delta t}
$$
$$
u(t) = K_p \cdot e(t) + K_i \cdot \text{积分项} + K_d \cdot \text{微分项}
$$

参数约束:

每个 PID 参数都被限制在指定的最小和最大边界内，以确保稳定性。

$$
Kp \in [Kp{\text{min}}, Kp_{\text{max}}]
$$
$$
Ki \in [Ki{\text{min}}, Ki_{\text{max}}]
$$
$$
Kd \in [Kd{\text{min}}, Kd_{\text{max}}]
$$

加热器系统动力学

模拟加热器的热行为，考虑热输入、环境损失和系统约束。

参数:

  • $K$: 加热器增益（热效率）。
  • $\tau$: 时间常数（响应速度）。
  • $T_d$: 死区时间（延迟步长）。
  • $T_s$: 采样时间。
  • $T_a$: 环境温度。
  • $h$: 热损失系数。
  • $P{\text{max}}$, $P{\text{min}}$: 功率约束。
  • $\Delta P_{\text{max}}$: 每步最大功率变化。

状态变量:

  • $y(t)$: 当前温度。
  • $u(t)$: 当前功率输入，受约束和延迟的影响。
  • $\text{delay_buffer}$: FIFO 缓冲区，用于模拟死区时间。

动态方程 (欧拉离散化):

$$
\frac{dy}{dt} = \frac{K \cdot u_{\text{delayed}}(t) - y(t)}{\tau} - h \cdot (y(t) - Ta)
$$
$$
y(t+1) = y(t) + \left( \frac{K \cdot u{\text{delayed}}(t) - y(t)}{\tau} - h \cdot (y(t) - T_a) \right) \cdot T_s
$$

控制信号处理:

  1. 限幅: 确保 $u(t)$ 位于 $[P{\text{min}}, P{\text{max}}]$ 内。
  2. 速率限制: 限制功率变化：

$$
u(t) = \text{clip}(u(t), u(t-1) - \Delta P{\text{max}}, u(t-1) + \Delta P{\text{max}})
$$

  1. 死区时间: 将 $u(t)$ 插入延迟缓冲区并检索 $u_{\text{delayed}}(t)$。

Smith预估器

Smith预估器通过使用无死区时间的系统模型来预测未来状态并相应地调整控制信号，从而预测系统延迟。

参数:

  • 与加热器系统共享相同的参数，除了 $T_d = 0$（无死区时间）。

预测方程:

$$
\frac{dy{\text{model}}}{dt} = \frac{K \cdot u(t) - y{\text{model}}(t)}{\tau} - h \cdot (y_{\text{model}}(t) - Ta)
$$
$$
y{\text{model}}(t+1) = y{\text{model}}(t) + \left( \frac{K \cdot u(t) - y{\text{model}}(t)}{\tau} - h \cdot (y_{\text{model}}(t) - T_a) \right) \cdot T_s
$$

补偿机制:

$$
u{\text{compensated}}(t) = u{\text{PID}}(t) + (SV(t) - y_{\text{model}}(t))
$$

其中 $SV(t)$ 是当前设定值，$y_{\text{model}}(t)$ 是预测温度。

仿真流程

仿真程序协调加热器系统、模糊PID控制器和Smith预估器在离散时间步长上的交互。

步骤:

  1. 初始化:

    • 设置初始 PID 参数 ($K{p0}, K{i0}, K_{d0}$)。
    • 定义设定值曲线，以固定间隔更改。
    • 初始化加热器系统、Smith 预估器、模糊 PID 控制器和 PID 控制器。
    • 初始化性能指标。

  2. 时间步进循环 ($t = 0$ 到 $t = \text{TIME_STEPS} - 1$):

    • 设定值分配: 分配当前设定值 $SV(t)$。
    • 过程变量 (PV): 检索当前温度 $y(t)$。
    • 误差计算:
      $$
      e(t) = SV(t) - PV(t)
      $$
      $$
      dE(t) = e(t) - e(t-1) \quad \text{(如果 } t > 0 \text{)}
      $$
    • 性能指标更新: 在设定值更改间隔，计算超调量、稳态误差和稳定时间。
    • 模糊化: 使用隶属度函数计算 $\muE$ 和 $\mu{dE}$。
    • 推理: 应用模糊规则确定 $\Delta K_p, \Delta K_i, \Delta K_d$。
    • 自适应 PID 调整: 更新 PID 参数：
      $$
      K_p(t) = \text{clip}(K_p(t-1) + \Delta Kp, Kp{\text{min}}, Kp_{\text{max}})
      $$
      $$
      K_i(t) = \text{clip}(K_i(t-1) + \Delta Ki, Ki{\text{min}}, Ki_{\text{max}})
      $$
      $$
      K_d(t) = \text{clip}(K_d(t-1) + \Delta Kd, Kd{\text{min}}, Kd_{\text{max}})
      $$
    • PID 控制输出:
      $$
      u_{\text{PID}}(t) = K_p(t) \cdot e(t) + K_i(t) \cdot \text{Integral}_e + K_d(t) \cdot dE(t)
      $$
    • Smith 预估器: 预测无延迟的 PV $\hat{y}(t)$。
    • 控制补偿:
      $$
      u{\text{compensated}}(t) = u{\text{PID}}(t) + (SV(t) - \hat{y}(t))
      $$
    • 加热器更新: 将 $u_{\text{compensated}}(t)$ 应用于加热器系统以获得 $y(t+1)$ 和 $u(t+1)$。
    • 数据记录: 存储 PV、SV、控制输出、PID 参数和性能指标。

  3. 最终性能指标: 循环结束后，计算最后一个设定值周期的指标。

性能指标

为了评估控制器的有效性，仿真跟踪每次设定值更改的以下指标：

  1. 超调量 (%):
    $$
    \text{超调量} = \left( \frac{\max(PV(t)) - SV}{SV} \right) \times 100
    $$
    测量 PV 超过设定值的程度。

  2. 稳态误差:
    $$
    \text{稳态误差} = \frac{1}{N} \sum_{t=T-N}^{T} (SV(t) - PV(t))
    $$
    响应最后 $N$ 步的平均误差。

  3. 稳定时间 (步长):
    PV 保持在设定值的 $\pm2\%$ 以内且之后不再超出此范围所需的步数。

用户界面

交互式图形界面，具有以下功能：

  1. 图表:

    • 系统响应: 显示 PV 和 SV 随时间的变化。
    • 控制信号: 显示控制输出 $u(t)$。
    • 加热器功率: 可视化应用于加热器的实际功率。
    • 自适应 PID 参数: 跟踪 $K_p, K_i, K_d$ 随时间的变化。
    • 性能指标: 绘制每次设定值更改的超调量、稳态误差和稳定时间。

  2. 滑块:

    • 初始 PID 参数:
      • Kp0: 初始比例增益 (0 到 50)。
      • Ki0: 初始积分增益 (0 到 10)。
      • Kd0: 初始微分增益 (0 到 10)。
    • 最大功率: 调整最大加热功率 (50 到 12,000)。

  3. 交互性:

    • 调整滑块会触发仿真程序使用更新的参数重新运行。
    • 图表动态更新以反映新的仿真结果。

结论

本文概述了与 Smith 预估器集成的模糊 PID 控制加热器系统仿真的算法结构和数学基础。该系统利用模糊逻辑进行自适应 PID 参数整定，在存在系统延迟和设定值变化的情况下增强控制性能。交互式仿真促进了控制动力学和性能指标的实时实验和可视化。

TODO：

  • 其他性能指标: 例如绝对误差积分 (IAE) 或平方误差积分 (ISE)。
  • 高级隶属度函数: 结合高斯或样条基隶属度函数以实现更平滑的过渡。
  • 优化算法: 使用优化算法来微调模糊规则库或 PID 参数范围。
  • 添加更多环境参数: 例如真空值、氩气/氟利昂等气体流量

附件

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import matplotlib.gridspec as gridspec

# 定义会员函数
def s_shaped(x, a, b):
    if x <= a:
        return 0
    elif a < x <= (a + b) / 2:
        return 2 * ((x - a) / (b - a)) ** 2
    elif (a + b) / 2 < x < b:
        return 1 - 2 * ((b - x) / (b - a)) ** 2
    else:
        return 1

def z_shaped(x, a, b):
    if x <= a:
        return 1
    elif a < x <= (a + b) / 2:
        return 1 - 2 * ((x - a) / (b - a)) ** 2
    elif (a + b) / 2 < x < b:
        return 2 * ((b - x) / (b - a)) ** 2
    else:
        return 0

def tri(x, a, b, c):
    if x <= a or x >= c:
        return 0
    elif a < x <= b:
        return (x - a) / (b - a)
    elif b < x < c:
        return (c - x) / (c - b)
    else:
        return 0

# 定义模糊控制器
class FuzzyPIDController:
    def __init__(self, E_max=100, dE_max=50):
        self.E_max = E_max
        self.dE_max = dE_max
        # 定义输出的独立映射数值
        self.output_map_Kp = {
            'NB': -0.5,
            'NM': -0.25,
            'NS': -0.1,
            'ZE': 0.0,
            'PS': 0.1,
            'PM': 0.25,
            'PB': 0.5
        }
        self.output_map_Ki = {
            'NB': -0.3,
            'NM': -0.15,
            'NS': -0.05,
            'ZE': 0.0,
            'PS': 0.05,
            'PM': 0.15,
            'PB': 0.3
        }
        self.output_map_Kd = {
            'NB': -0.2,
            'NM': -0.1,
            'NS': -0.05,
            'ZE': 0.0,
            'PS': 0.05,
            'PM': 0.1,
            'PB': 0.2
        }
        # 定义规则表
        self.rules = [
            # NB row
            [('PB', 'ZE', 'NS'), ('PM', 'PS', 'ZE'), ('PM', 'PS', 'PS'),
             ('PM', 'PM', 'PS'), ('PS', 'PM', 'PS'), ('ZE', 'PM', 'PM'),
             ('ZE', 'PS', 'PM')],
            # NM row
            [('PM', 'PS', 'NS'), ('PM', 'PS', 'PS'), ('PM', 'PM', 'PS'),
             ('PS', 'PM', 'PS'), ('ZE', 'PM', 'PM'), ('ZE', 'PS', 'PM'),
             ('NS', 'ZE', 'PM')],
            # NS row
            [('PM', 'PS', 'PS'), ('PM', 'PM', 'PS'), ('PS', 'PM', 'PS'),
             ('ZE', 'PM', 'PM'), ('ZE', 'PS', 'PM'), ('NS', 'ZE', 'PM'),
             ('NM', 'NS', 'PB')],
            # ZE row
            [('PM', 'PM', 'PS'), ('PS', 'PM', 'PS'), ('ZE', 'PM', 'PM'),
             ('ZE', 'ZE', 'ZE'), ('ZE', 'NS', 'PM'), ('NM', 'NS', 'PB'),
             ('NM', 'NM', 'PB')],
            # PS row
            [('PS', 'PM', 'PS'), ('ZE', 'PM', 'PM'), ('ZE', 'PS', 'PM'),
             ('NS', 'ZE', 'PM'), ('NM', 'NS', 'PB'), ('NM', 'NM', 'PB'),
             ('NB', 'NB', 'PB')],
            # PM row
            [('ZE', 'PM', 'PM'), ('ZE', 'PS', 'PM'), ('NS', 'ZE', 'PM'),
             ('NM', 'NS', 'PB'), ('NM', 'NM', 'PB'), ('NB', 'NB', 'PB'),
             ('NB', 'NB', 'PB')],
            # PB row
            [('ZE', 'PS', 'PM'), ('NS', 'ZE', 'PM'), ('NM', 'NS', 'PB'),
             ('NM', 'NM', 'PB'), ('NB', 'NB', 'PB'), ('NB', 'NB', 'PB'),
             ('NB', 'NB', 'PB')]
        ]
        self.E_labels = ['NB', 'NM', 'NS', 'ZE', 'PS', 'PM', 'PB']
        self.dE_labels = ['NB', 'NM', 'NS', 'ZE', 'PS', 'PM', 'PB']

    def fuzzify(self, E, dE):
        mu_E = {}
        mu_dE = {}
        # Fuzzify E
        mu_E['NB'] = z_shaped(E, -self.E_max, -self.E_max / 2)
        mu_E['NM'] = tri(E, -self.E_max, -self.E_max / 2, 0)
        mu_E['NS'] = tri(E, -self.E_max / 2, 0, self.E_max / 4)
        mu_E['ZE'] = tri(E, -self.E_max / 4, 0, self.E_max / 4)
        mu_E['PS'] = tri(E, 0, self.E_max / 4, self.E_max / 2)
        mu_E['PM'] = tri(E, self.E_max / 3, self.E_max / 2, 2 * self.E_max / 3)
        mu_E['PB'] = s_shaped(E, self.E_max / 2, self.E_max)
        # Fuzzify dE
        mu_dE['NB'] = z_shaped(dE, -self.dE_max, -self.dE_max / 2)
        mu_dE['NM'] = tri(dE, -self.dE_max, -self.dE_max / 2, 0)
        mu_dE['NS'] = tri(dE, -self.dE_max / 2, 0, self.dE_max / 2)
        mu_dE['ZE'] = tri(dE, -self.dE_max / 4, 0, self.dE_max / 4)
        mu_dE['PS'] = tri(dE, 0, self.dE_max / 4, self.dE_max / 2)
        mu_dE['PM'] = tri(dE, self.dE_max / 3, self.dE_max / 2, 2 * self.dE_max / 3)
        mu_dE['PB'] = s_shaped(dE, self.dE_max / 2, self.dE_max)
        return mu_E, mu_dE

    def infer(self, mu_E, mu_dE):
        # 初始化输出会员度列表
        delta_Kp = []
        delta_Ki = []
        delta_Kd = []

        # 遍历所有规则
        for i, E_label in enumerate(self.E_labels):
            for j, dE_label in enumerate(self.dE_labels):
                alpha = min(mu_E[E_label], mu_dE[dE_label])  # 取最小
                if alpha > 0:
                    rule = self.rules[i][j]
                    delta_Kp.append((self.output_map_Kp[rule[0]], alpha))
                    delta_Ki.append((self.output_map_Ki[rule[1]], alpha))
                    delta_Kd.append((self.output_map_Kd[rule[2]], alpha))

        # 聚合: 对每个输出变量取加权平均
        def weighted_average(output):
            numerator = sum([val * alpha for val, alpha in output])
            denominator = sum([alpha for _, alpha in output])
            if denominator == 0:
                return 0
            return numerator / denominator

        delta_Kp_value = weighted_average(delta_Kp)
        delta_Ki_value = weighted_average(delta_Ki)
        delta_Kd_value = weighted_average(delta_Kd)

        return delta_Kp_value, delta_Ki_value, delta_Kd_value

# 定义PID控制器
class PIDController:
    def __init__(self, Kp=2.0, Ki=0.5, Kd=0.1, 
                 Kp_min=0, Kp_max=10, 
                 Ki_min=0, Ki_max=10, 
                 Kd_min=0, Kd_max=10):
        self.Kp = Kp
        self.Ki = Ki
        self.Kd = Kd
        self.Kp_min = Kp_min
        self.Kp_max = Kp_max
        self.Ki_min = Ki_min
        self.Ki_max = Ki_max
        self.Kd_min = Kd_min
        self.Kd_max = Kd_max
        self.integral = 0
        self.prev_error = 0

    def update_parameters(self, delta_Kp, delta_Ki, delta_Kd):
        self.Kp = np.clip(self.Kp + delta_Kp, self.Kp_min, self.Kp_max)
        self.Ki = np.clip(self.Ki + delta_Ki, self.Ki_min, self.Ki_max)
        self.Kd = np.clip(self.Kd + delta_Kd, self.Kd_min, self.Kd_max)

    def compute(self, error, dt=1):
        # 积分
        self.integral += error * dt
        # 微分
        derivative = (error - self.prev_error) / dt
        self.prev_error = error

        return self.Kp * error + self.Ki * self.integral + self.Kd * derivative

# 定义加热器系统模型，包含环境热量散失、加热功率限制和加热功率变化率限制
class HeaterSystem:
    def __init__(self, K=1.0, tau=10.0, Td=5, Ts=0.05, ambient_temperature=20.0, heat_loss_coefficient=0.1, 
                 max_power=100.0, min_power=0.0, power_rate_limit=10.0):
        self.K = K
        self.tau = tau
        self.Td = Td  # 滞后时间（延迟步数）
        self.Ts = Ts    # 采样时间
        self.ambient_temperature = ambient_temperature  # 环境温度
        self.heat_loss_coefficient = heat_loss_coefficient  # 热量散失系数
        self.delay_buffer = [0] * Td  # 滞后缓冲区
        self.y = ambient_temperature  # 初始温度
        self.max_power = max_power  # 最大加热功率
        self.min_power = min_power  # 最小加热功率
        self.power_rate_limit = power_rate_limit  # 加热功率变化率限制
        self.prev_delayed_u = 0.0  # 记录上一个延迟的加热功率

    def update(self, u):
        # 限制控制信号在加热功率范围内
        u = np.clip(u, self.min_power, self.max_power)

        # 限制加热功率变化率
        delta_u = u - self.prev_delayed_u
        delta_u = np.clip(delta_u, -self.power_rate_limit, self.power_rate_limit)
        u = self.prev_delayed_u + delta_u

        # 添加控制输入到延迟缓冲区
        self.delay_buffer.append(u)
        delayed_u = self.delay_buffer.pop(0)
        self.prev_delayed_u = delayed_u  # 更新上一个延迟的加热功率

        # 欧拉法离散化一阶惯性环节，加上环境热量散失
        dy = (self.K * delayed_u - self.y) * self.Ts / self.tau
        # 环境热量散失
        dy -= (self.y - self.ambient_temperature) * self.heat_loss_coefficient * self.Ts
        self.y += dy

        return self.y, delayed_u  # 返回延迟后的加热功率

# 定义史密斯预测器
class SmithPredictor:
    def __init__(self, model: HeaterSystem, Td=5):
        self.model = model
        self.Td = Td
        self.model_delay_buffer = [0] * Td
        self.model_output = self.model.ambient_temperature  # 初始预测输出

    def predict(self, u):
        # 限制控制信号在加热功率范围内
        u = np.clip(u, self.model.min_power, self.model.max_power)

        # 预测模型不含滞后
        dy = (self.model.K * u - self.model_output) * self.model.Ts / self.model.tau
        # 环境热量散失
        dy -= (self.model_output - self.model.ambient_temperature) * self.model.heat_loss_coefficient * self.model.Ts
        self.model_output += dy
        return self.model_output

# 定义仿真函数
def simulate(Kp0, Ki0, Kd0, E_max=100, dE_max=50, TIME_STEPS=2000, 
             SETPOINT_CHANGE_INTERVAL=500, max_power=12000.0, 
             min_power=0.0, power_rate_limit=100.0, heater_power=0):
    # 设定点配置：每SETPOINT_CHANGE_INTERVAL步改变一次设定点
    SETPOINT_VALUES = [200, 400, 300, 600]
    num_changes = TIME_STEPS // SETPOINT_CHANGE_INTERVAL
    SETPOINT_PROFILE = np.concatenate([np.ones(SETPOINT_CHANGE_INTERVAL) * SETPOINT_VALUES[i % len(SETPOINT_VALUES)] for i in range(num_changes)])
    if len(SETPOINT_PROFILE) < TIME_STEPS:
        SETPOINT_PROFILE = np.concatenate([SETPOINT_PROFILE, np.ones(TIME_STEPS - len(SETPOINT_PROFILE)) * SETPOINT_VALUES[-1]])

    # 初始化模块
    heater = HeaterSystem(K=1, tau=10.0, Td=1, Ts=0.1, ambient_temperature=20.0, 
                          heat_loss_coefficient=0.1, max_power=max_power, min_power=min_power, 
                          power_rate_limit=power_rate_limit)
    smith_model = HeaterSystem(K=1, tau=10.0, Td=0, Ts=0.1, ambient_temperature=20.0, 
                               heat_loss_coefficient=0.1, max_power=max_power, min_power=min_power, 
                               power_rate_limit=power_rate_limit)  # 模型不含滞后
    smith = SmithPredictor(model=smith_model, Td=1)
    fuzzy_controller = FuzzyPIDController(E_max=E_max, dE_max=dE_max)
    pid = PIDController(Kp=Kp0, Ki=Ki0, Kd=Kd0, 
                        Kp_min=0, Kp_max=50, 
                        Ki_min=0, Ki_max=10, 
                        Kd_min=0, Kd_max=10)

    # 初始化记录变量
    PV_history = []
    SV_history = []
    error_history = []
    control_history = []
    heater_power_history = []
    Kp_history = []
    Ki_history = []
    Kd_history = []
    overshoot_history = []
    steady_state_error_history = []
    settling_time_history = []

    # 初始化性能指标
    current_setpoint = SETPOINT_PROFILE[0]
    response_start_index = 0
    overshoot = 0
    steady_state_error = 0
    settling_time = 0
    settling_time_threshold = 0.02  # 2%
    within_settle = False

    # 仿真循环
    for t in range(TIME_STEPS):
        SV = SETPOINT_PROFILE[t]
        PV = heater.y
        error = SV - PV
        dE = error - (error_history[-1] if t > 0 else 0)

        # 检测设定点变化
        if t % SETPOINT_CHANGE_INTERVAL == 0 and t != 0:
            # 计算并记录性能指标
            # 超调量
            max_PV = max(PV_history[response_start_index:t]) if t > response_start_index else PV
            overshoot = ((max_PV - current_setpoint) / current_setpoint) * 100 if current_setpoint != 0 else 0
            overshoot_history.append(overshoot)

            # 稳态误差 (取最后50步的平均误差)
            steady_state_error = np.mean(error_history[max(t-50, response_start_index):t]) if t > response_start_index else error
            steady_state_error_history.append(steady_state_error)

            # 调节时间 (找到第一次进入并保持在±2%范围内的时间步)
            target_min = current_setpoint * (1 - settling_time_threshold)
            target_max = current_setpoint * (1 + settling_time_threshold)
            settling_time = 0
            for i in range(response_start_index, t):
                if target_min <= PV_history[i] <= target_max:
                    # Check if it stays within the range for the remaining steps
                    if all(target_min <= pv <= target_max for pv in PV_history[i:t]):
                        settling_time = i - response_start_index
                        break
            settling_time_history.append(settling_time)

            # 更新当前设定点和响应起始索引
            current_setpoint = SV
            response_start_index = t

        # 模糊控制器调整PID参数
        mu_E, mu_dE = fuzzy_controller.fuzzify(error, dE)
        delta_Kp, delta_Ki, delta_Kd = fuzzy_controller.infer(mu_E, mu_dE)

        if heater_power != 0:
            delta_Kp *= (heater_power / max_power) # 根据加热功率调整PID参数

        pid.update_parameters(delta_Kp, delta_Ki, delta_Kd)

        # PID计算控制输出
        u_pid = pid.compute(error, dt=heater.Ts)

        # 史密斯预测器预测PV
        predicted_PV = smith.predict(u_pid)

        # 控制输出补偿
        u_compensated = u_pid + (SV - predicted_PV)

        # 更新加热器系统，并获取延迟后的加热功率
        y, heater_power = heater.update(u_compensated)

        # 记录数据
        PV_history.append(y)
        SV_history.append(SV)
        error_history.append(error)
        control_history.append(u_compensated)
        heater_power_history.append(heater_power)
        Kp_history.append(pid.Kp)
        Ki_history.append(pid.Ki)
        Kd_history.append(pid.Kd)

    # 最后一个设定点的性能指标
    max_PV = max(PV_history[response_start_index:]) if TIME_STEPS > response_start_index else PV_history[-1]
    overshoot = ((max_PV - current_setpoint) / current_setpoint) * 100 if current_setpoint != 0 else 0
    overshoot_history.append(overshoot)
    steady_state_error = np.mean(error_history[response_start_index:]) if TIME_STEPS > response_start_index else error_history[-1]
    steady_state_error_history.append(steady_state_error)
    settling_time = 0
    target_min = current_setpoint * (1 - settling_time_threshold)
    target_max = current_setpoint * (1 + settling_time_threshold)
    for i in range(response_start_index, TIME_STEPS):
        if target_min <= PV_history[i] <= target_max:
            if all(target_min <= pv <= target_max for pv in PV_history[i:]):
                settling_time = i - response_start_index
                break
    settling_time_history.append(settling_time)

    return (PV_history, SV_history, control_history, heater_power_history, 
            Kp_history, Ki_history, Kd_history, 
            overshoot_history, steady_state_error_history, settling_time_history)

# 初始化仿真参数
INITIAL_Kp = 2.0
INITIAL_Ki = 0.5
INITIAL_Kd = 0.1

# 运行初始仿真
SETPOINT_CHANGE_INTERVAL=500
results = simulate(INITIAL_Kp, INITIAL_Ki, INITIAL_Kd, TIME_STEPS=2000, 
                   SETPOINT_CHANGE_INTERVAL=SETPOINT_CHANGE_INTERVAL, 
                   max_power=12000, min_power=0, power_rate_limit=10)
PV_history, SV_history, control_history, heater_power_history, Kp_history, Ki_history, Kd_history, \
overshoot_history, steady_state_error_history, settling_time_history = results

# 创建图形和子图，增加一个子图用于展示加热功率和性能指标
fig = plt.figure(figsize=(15, 20))
gs = gridspec.GridSpec(5, 1, height_ratios=[2, 1, 1, 1, 1], hspace=0.4)

# 子图1：系统响应
ax1 = plt.subplot(gs[0])
line1_PV, = ax1.plot(PV_history, label='Process Variable (PV)', color='b')
line1_SV, = ax1.plot(SV_history, label='Setpoint (SV)', linestyle='--', color='r')
ax1.set_ylabel('Temperature')
ax1.set_title('System Response')
ax1.legend()
ax1.grid(True)

# 子图2：控制输出
ax2 = plt.subplot(gs[1], sharex=ax1)
line2, = ax2.plot(control_history, label='Control Output (u)', color='g')
ax2.set_ylabel('Control Output')
ax2.set_title('Control Signal')
ax2.legend()
ax2.grid(True)

# 子图3：加热功率
ax3 = plt.subplot(gs[2], sharex=ax1)
line3, = ax3.plot(heater_power_history, label='Heater Power (u)', color='orange')
ax3.set_ylabel('Heater Power')
ax3.set_title('Heater Power')
ax3.legend()
ax3.grid(True)

# 子图4：PID参数
ax4 = plt.subplot(gs[3], sharex=ax1)
line4_Kp, = ax4.plot(Kp_history, label='Kp', color='m')
line4_Ki, = ax4.plot(Ki_history, label='Ki', color='c')
line4_Kd, = ax4.plot(Kd_history, label='Kd', color='y')
ax4.set_ylabel('PID Parameters')
ax4.set_title('Adaptive PID Parameters')
ax4.legend()
ax4.grid(True)

# 子图5：性能指标
ax5 = plt.subplot(gs[4])
steps = [SETPOINT_CHANGE_INTERVAL * i for i in range(len(overshoot_history))]
ax5.plot(steps, overshoot_history, label='Overshoot (%)', marker='o', color='purple')
ax5.plot(steps, steady_state_error_history, label='Steady-State Error', marker='x', color='brown')
ax5.plot(steps, settling_time_history, label='Settling Time (steps)', marker='^', color='black')
ax5.set_xlabel('Time Step')
ax5.set_ylabel('Performance Metrics')
ax5.set_title('Performance Metrics per Setpoint Change')
ax5.legend()
ax5.grid(True)

# 调整子图间距以留出滑块空间
plt.subplots_adjust(bottom=0.35)

# 添加滑块
axcolor = 'lightgoldenrodyellow'
ax_Kp = plt.axes([0.15, 0.25, 0.65, 0.03], facecolor=axcolor)
ax_Ki = plt.axes([0.15, 0.18, 0.65, 0.03], facecolor=axcolor)
ax_Kd = plt.axes([0.15, 0.11, 0.65, 0.03], facecolor=axcolor)
ax_max_power = plt.axes([0.15, 0.04, 0.65, 0.03], facecolor=axcolor)

slider_Kp = Slider(ax_Kp, 'Kp0', 0.0, 50.0, valinit=INITIAL_Kp, valstep=0.1)
slider_Ki = Slider(ax_Ki, 'Ki0', 0.0, 10.0, valinit=INITIAL_Ki, valstep=0.1)
slider_Kd = Slider(ax_Kd, 'Kd0', 0.0, 10.0, valinit=INITIAL_Kd, valstep=0.1)
slider_max_power = Slider(ax_max_power, 'Max Power', 50.0, 12000.0, valinit=2000.0, valstep=1.0)

# 定义更新函数
def update(val):
    Kp0 = slider_Kp.val
    Ki0 = slider_Ki.val
    Kd0 = slider_Kd.val
    max_power = slider_max_power.val
    min_power = 0.0  # 可以根据需要调整最小加热功率
    power_rate_limit = 500.0  # 固定变化率限制，可扩展为可调

    # 重新运行仿真
    results = simulate(Kp0, Ki0, Kd0, TIME_STEPS=2000, 
                      SETPOINT_CHANGE_INTERVAL=SETPOINT_CHANGE_INTERVAL, 
                      max_power=max_power, min_power=min_power, 
                      power_rate_limit=power_rate_limit)
    PV, SV, u, heater_power, Kp, Ki, Kd, overshoot, steady_state_error, settling_time = results

    # 更新系统响应图
    line1_PV.set_ydata(PV)
    line1_SV.set_ydata(SV)
    ax1.relim()
    ax1.autoscale_view()

    # 更新控制输出图
    line2.set_ydata(u)
    ax2.relim()
    ax2.autoscale_view()

    # 更新加热功率图
    line3.set_ydata(heater_power)
    ax3.relim()
    ax3.autoscale_view()

    # 更新PID参数图
    line4_Kp.set_ydata(Kp)
    line4_Ki.set_ydata(Ki)
    line4_Kd.set_ydata(Kd)
    ax4.relim()
    ax4.autoscale_view()

    # 更新性能指标图
    steps = [500 * i for i in range(len(overshoot))]
    ax5.clear()
    ax5.plot(steps, overshoot, label='Overshoot (%)', marker='o', color='purple')
    ax5.plot(steps, steady_state_error, label='Steady-State Error', marker='x', color='brown')
    ax5.plot(steps, settling_time, label='Settling Time (steps)', marker='^', color='black')
    ax5.set_xlabel('Time Step')
    ax5.set_ylabel('Performance Metrics')
    ax5.set_title('Performance Metrics per Setpoint Change')
    ax5.legend()
    ax5.grid(True)

    fig.canvas.draw_idle()

# 注册滑块更新事件
slider_Kp.on_changed(update)
slider_Ki.on_changed(update)
slider_Kd.on_changed(update)
slider_max_power.on_changed(update)

plt.show()
此作者没有提供个人介绍
最后更新于 2024-10-14
啥也没有呀