2022年卡尔曼滤波器在热力系统中的应用.docx

精选学习资料 - - - - - - - - - 5 热 力 系 统 的 设 计本 文 应 用 卡 尔 曼 滤 波 对 温 室 室 内 温 度x 1t和室 外 温 度x2t,进 行 估 计 ,我们用 测 温 传 感 器 对 室 内 温 度x 1t进行测量,并采 用A/D转 换 器 对 测 量 值 进 行等 时 间 间 隔 取 样 , 取 样 时 间 间 隔 为 T ,时 标 为 k ；热力系统如下示wk温 室 热 力 系 统A/D vky kx2tx 1t图 5-1 5.1 建 立 数 学 模 型上 述 系 统 的 离 散 状 态 模 型 如 下 图ukBkw kXk1Z1XkCkv k kyZA k图 5-2 离 散 系 统 状 态 空 间 模 型此 热 力 系 统 的 状 态 方 程 为 ：名师归纳总结 x 1k1.06270.361x 1k.00125w k（5-1）第 1 页,共 12 页x 2k1 .00900. 833x 2k.00575- - - - - - -精选学习资料 - - - - - - - - - 量 测 方 程 为 ：yk1 w10 x 1k1 v k1（5-2）是在x2k1上 二 式 中 ：k是 因 环 境 温 度 随 机 变 化 而 产 生 的 过 程 噪 声 ；vk1 测 量 过 程 中 引 入 的 量 测 噪 声 ；当 k 变 化 时 ,他 们 将 构 成 两 个 互 不 相 关 的 白 噪 声 ；当 k 变 化 时 , 它 们 将 构 成 两 个 互 不 相 关 的 白 噪 声 序 列 wk和 vk,其 方 差分 别 为 ：Ewk209Evk2. 2我 们 可 根 据 此 求 出 过 程 噪 声 和 量 测 噪 声 的 协 方 差 矩 阵QE.00125 w k220. 01250. 0575Ewk20. 05750 .0125 EwkE0 .0575 wk2R.0 25.2 状 态 观 测 器 的 设 计由 于 本 文 中 的 系 统 方 程 和 量 测 方 程 的 形 式 是xk1 Axkwk1y k1Cx k1v k而 不 是xkAxk1 kwk1ykCx kv故 为 使 时 标 k 的取法一样,卡尔 曼 滤 波 器 的 递 推 公 式 可 写 成 如 下名师归纳总结 x k1 Ax kKk1 yk1 CAxkk15-3 第 2 页,共 12 页5-4 Kk1 P 1 k1CTCP 1k1 CTR 1 P 1k1 kT Ak15-5 APQP k1P 1 k1Kk1 CP 1 k5-6 - - - - - - -精选学习资料 - - - - - - - - - wkvk1 ACyk1+ 1CA+ K k1+ 1 CAx kxk图 5-3 状 态 观 测 图5.3 卡 尔 曼 滤 波 器 的 初 始 值 叠 代如 假 定 我 们 已 经 知 道 k =0 时 内 外 室 的 初 始 温 度就1P00,K0 0,Px 0 20250 0利 用 （ 5-5 ） 式 , 可 得P 1 1 AP0ATQQ.0 0 0 1 4 .0 0 0 6 5.0 0 0 6 5 .0 0 2 9 8利 用 （ 5-4 ） 式 , 可 得K1 P 1 1 CTCP 1 1CTR1100 .00140 . 006510 . 210 . 00140 . 006510 . 00650 . 02980.000650 . 02980.00069.0 0323再 利 用 （ 5-6 ） 式 , 可 得名师归纳总结 P 1 P 11K1 CP 1 1第 3 页,共 12 页- - - - - - -精选学习资料 - - - - - - - - - 0 . 00140 . 0065.00069100 . 00140 . 00650 . 00650 . 02980. 0323.0 00650 . 02980 . 0 0 1 4 0 . 0 0 6 4、0 . 0 0 6 4 0 . 0 2 9 8如 将 上 述 递 推 过 程 继 续 进 行 下 去 , 我 们 即 可 依 次 求 出 k =110 时P 1kKk、Pk值；k =1 10 时 的 值 卡 尔 曼 滤 波 增 益 值 如 图 5-9示 ；从 图 5-9中 可 知 , 当 卡 尔 曼 滤 波 增 益 在 k=8 时 已 基 本 达 到 其 稳 定 值k0 . 0190 . 018可 写 出 稳 态 时 的 卡 尔 曼 滤 波 器 的 状 态 估 计 方 程 ：x 1k1 0 .6270. 361x 1k00.6270.361x 1k（ 5-7）x2k1 0 .0900.833x2k0. 200yk1 10.2480. 0910. 833x2k或写成x 1 k 1 0 . 627 0 . 361 x 1 k x 2 k 1 0 . 090 0 . 833 x 2 k 0 . 200y k 1 0 . 627 x 1 k 0 . 361 x 2 k （5-8）0 . 248显 然 , 只 要 选 定 对 系 统 状 态 的 初 始 估 计 值 x 0 , 即 可 利 用 依 次 求 得 的 滤波 增 益 值 K k 和 滤 波 估 计 方 程（ 5-3 ）式 对 系 统 状 态 进 行 递 推 估 计 ；但 是 ,由于 在 k=10 后 K k 值 基 本 趋 于 稳 定 ,故 也 可 利 用 稳 态 滤 波 估 计 方 程（ 5-8 ）式 来对 系 统 状 态 进 行 递 推 估 计 ； 此 时 , 由 于 在 最 初 几 次 测 量 中 采 用 了 较 高 的 滤 波增 益 稳 态 值 , 对 系 统 状 态 的 前 几 次 估 计 会 出 现 较 大 的 误 差 ； 当 然 , 大 约 在 十个 取 样 周 期 之 后 , 这 种 情 况 即 可 消 除 ；名师归纳总结 将 上 述 结 果 通 过 MATLAB运 行 之 后 的 图 如 图 5.5-5.8所 示程 序 见 附 录 ：第 4 页,共 12 页- - - - - - -精选学习资料 - - - - - - - - - 25Figure 1 - Temperature True, Measured, and Estimated20c perature Tem151050-5012345678910Time sec图 5-5 真 实 值 、 测 量 值 及 预 测 估 计 值 随 时 间 变 化 的 图在 图 5-5 中 ：红 色 -温 度 的 真 实 值 ；蓝 色 -温 度 的 预 测 估 计 值 ；绿 色 -温 度 的 测 量 值 ；由 图 5-5 可 看 出 真 实 值 与 预 测 估 计 值 变 化 基 本 一 致 , 可 看 成 近 似 相 等 ,而 测 量 值 因 受 到 量 测 噪 声 影 响 而 有 所 波 动 变 化 , 不 过 取 值 在 真 实 值 的 小 范 围内 上 下 波 动 , 对 以 后 的 预 测 估 计 影 响 较 小 ；名师归纳总结 - - - - - - -第 5 页,共 12 页精选学习资料 - - - - - - - - - Figure 2 - Temperature Measurement Error and Temperature Estimation Error 1.51Errorc perature Tem0.50-0.5-1012345678910Time sec图 5-6将 图 5-5图 5-6 测 量 值 与 预 测 估 计 值 的 误 差中 的 预 测 估 计 值 与 测 量 值 的 波 动 曲 线 进 行 了 放 大 , 由 此 ,我 们 可 以 看 出 , 预 测 估 计 值 的 误 差 基 本 上 趋 向 于 零 , 说 明 预 测 估 计 值 与 真 实值 趋 于 相 等 , 可 见 卡 尔 曼 滤 波 器 预 测 估 计 的 准 确 性 ；而 测 量 值 误 差 由 于 受 到 量 测 噪 声 的 影 响 就 在 -11 化 , 波 动 幅 度 较 小 , 对 测 量 结 果 影 响 不 大 ；摄 氏 度 之 间 上 下 波 动 变名师归纳总结 - - - - - - -第 6 页,共 12 页精选学习资料 - - - - - - - - - 25Figure 3 - Temperature True and Estimated20c/sec perature Tem151050-5012345678910Time sec图 5-7 真实值与猜测估量值的变化曲线图 5-7 中说明白真实值与猜测估量值的变化趋于一样,近似相等；0.1Figure 4 -Temperature Estimation Error0.05Errorc/sec perature Tem0-0.05-0.1-0.15012345678910Time sec名师归纳总结 图 5-8 温度猜测估量值误差的波动曲线第 7 页,共 12 页- - - - - - -精选学习资料 - - - - - - - - - 从图 5-8 中我们可以看到猜测估量值误差的详细波动情形,它的波动变化始终保持在-0.120.1 摄氏度之间,从温度猜测估量的角度来说,此误差对整体的猜测估量值基本没有影响,从而使得猜测估量值与真实值时间的误差基本趋于零；图 5-9 卡 尔 曼 滤 波 增 益 变 化 图由 图 5-9可 知 , 当 k=8 时 , 增 益 Kk逐 渐 趋 向 于 稳 定 ； 蓝 线 是 温 度x 1t的 增益 , 绿 线 是 温 度x 2t的 增 益附 录function kalmanduration, dt % function kalmanduration, dt % Kalman filter simulation for a vehicle travelling along a road. % INPUTS 名师归纳总结 - - - - - - -第 8 页,共 12 页精选学习资料 - - - - - - - - - % duration = length of simulation seconds % dt = step size seconds measnoise = 10; % position measurement noise feet accelnoise = 0.2; % acceleration noise feet/sec2 a = 0.627 0.361;0.090 0.833; % transition matrix b = 0;0; % input matrix c = 1 0; % measurement matrix x = 20;25; % initial state vector xhat = x; % initial state estimate Sz = measnoise2; % measurement error covariance Sw = accelnoise2 * dt4/4 dt3/2; dt3/2 dt2; % process noise cov P = Sw; % initial estimation covariance % Initialize arrays for later plotting. pos = ; % true position array poshat = ; % estimated position array posmeas = ; % measured position array vel = ; % true velocity array velhat = ; % estimated velocity array for t = 0 : dt: duration, % Use a constant commanded acceleration of 1 foot/sec2. 名师归纳总结 - - - - - - -第 9 页,共 12 页精选学习资料 - - - - - - - - - u = 0; % Simulate the linear system. ProcessNoise = accelnoise * dt2/2*randn; dt*randn; x = a * x + b * u + ProcessNoise; % Simulate the noisy measurement MeasNoise = measnoise * randn; y = c * x + MeasNoise; % Extrapolate the most recent state estimate to the present time. xhat = a * xhat + b * u; % Form the Innovation vector. Inn = y - c * xhat; % Compute the covariance of the Innovation. s = c * P * c + Sz;% Form the Kalman Gain matrix. K = a * P * c * invs;% Update the state estimate. xhat = xhat + K * Inn; % Compute the covariance of the estimation error. P = a * P * a - a * P * c * invs * c * P * a + Sw;% Save some parameters for plotting later. pos = pos; x1; 名师归纳总结 - - - - - - -第 10 页,共 12 页精选学习资料 - - - - - - - - - posmeas = posmeas; y; poshat = poshat; xhat1; vel = vel; x2; velhat = velhat; xhat2; end % Plot the results close all; t = 0 : dt : duration; figure; plott,pos, t,posmeas, t,poshat; grid; xlabel Time sec ;ylabel Position feet;titleFigure 1 - Vehicle Position True, Measured, and Estimated figure; plott,pos-posmeas, t,pos-poshat; grid; 名师归纳总结 xlabel Time sec ;第 11 页,共 12 页ylabel Position Error feet;titleFigure 2 - Position Measurement Error and Position Estimation Error;- - - - - - -精选学习资料 - - - - - - - - - figure; plott,vel, t,velhat; grid; xlabel Time sec ;ylabel Velocity feet/sec;titleFigure 3 - Velocity True and Estimatedfigure; plott,vel-velhat; grid; 名师归纳总结 xlabel Time sec ;第 12 页,共 12 页ylabel Velocity Error feet/sectitleFigure 4 - Velocity Estimation Error- - - - - - -