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

《2022年卡尔曼滤波器在热力系统中的应用 .pdf》由会员分享，可在线阅读，更多相关《2022年卡尔曼滤波器在热力系统中的应用 .pdf（12页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、5 热 力 系 统 的 设 计本 文 应 用 卡 尔 曼 滤 波 对 温 室 室 内 温 度)(1tx和室 外 温 度)(2tx，进 行 估 计 ，我们用 测 温 传 感 器 对 室 内 温 度)(1tx进行测量，并采 用A/D转 换 器 对 测 量 值 进 行等 时 间 间 隔 取 样 ， 取 样 时 间 间 隔 为T，时 标 为k。热力系统如下示图 5-1 温 室 热 力 系 统5.1 建 立 数 学 模型上 述 系 统 的 离 散 状 态 模 型 如 下 图图 5-2 离 散 系 统 状 态 空 间 模 型此 热 力 系 统 的 状 态 方 程 为 ：)(0125. 00575. 0)(

2、)(0.361627. 0833.0090. 0)1()1(1212kwkxkxkxkx（5-1）)(1txA/D )(ky)(kv)(kw)(2txkB1ZZkCkAkX1kXku)(kvkykw精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 12 页量 测 方 程 为 ：) 1(01)1()1()1(12kvkykxkx（5-2）上 二 式 中 ：)(kw是 因 环 境 温 度 随 机 变 化 而 产 生 的 过 程 噪 声 ；)1(kv是在测 量 过 程 中 引 入 的 量 测 噪 声 。 当k变 化 时 ， 他 们 将 构 成

3、两 个 互 不 相 关 的 白 噪 声 。当k变 化 时 ， 它 们 将 构 成 两 个 互 不 相 关 的 白 噪 声 序 列 )(kw和 )(kv，其 方 差分 别 为 ：9)(2kwE2.0)(2kvE我 们 可 根 据 此 求 出 过 程 噪 声 和 量 测 噪 声 的 协 方 差 矩 阵2222)(0575.0)(0125.00575.0)(0575.00125.0)(0125. 0kwEkwEkwEkwEQ2. 0R5.2 状 态 观 测 器的 设 计由 于 本 文 中 的 系 统 方 程 和 量 测 方 程 的 形 式 是)()()1(kwkAxkx) 1() 1() 1(kvk

4、Cxky而 不 是) 1()1()(kwkAxkx)()()(kvkCxky故 为 使 时 标k的取法一致，卡尔 曼 滤 波 器 的 递 推 公 式 可 写 成 如 下)()1()1()()1(kxCAkykKkxAkx(5-3) 111)1()1() 1()1(kRCkCPCkPkKTT(5-4) )()()1(1kQAkAPkPT(5-5) ) 1()1() 1() 1(11kCPkKkPkP(5-6) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 12 页图 5-3 状 态 观 测 图5.3 卡 尔 曼 滤 波器 的 初 始 值

5、 叠 代若 假 定 我 们 已 经 知 道k=0 时 内 外 室 的 初 始 温 度2520)0(x则0)0(1P,0)0(K,0)0(P利 用 （ 5-5 ） 式 ， 可 得QAAPPT)0()1(1Q0 2 98. 00 0 6 5. 00 0 65. 00 0 1 4. 0利 用 （ 5-4 ） 式 ， 可 得111) 1()1()1 (RCCPCPKTT12.0010298.00065. 00065.00014.001010298.00065.00065.00014.00323. 00069. 0再 利 用 （ 5-6 ） 式 ， 可 得) 1()1 () 1 ()1(11CPKPPC

6、A1CA) 1(kKAC)(kw)1(kv)(kx+ + CA+ ) 1(ky)1(kx精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 12 页0298.00065. 00065.00014.0010323.00069. 00298.00065.00065.00014.00 2 9 8.00 0 6 4.00 0 6 4.00 0 1 4.0若 将 上 述 递 推 过 程 继 续 进 行 下 去 ， 我 们 即 可 依 次 求 出k=110 时)(1kP、)(kK、)(kP值。k=1 10 时 的 值 卡 尔 曼 滤 波 增 益 值 如

7、 图 5-9示 。从 图 5-9中 可 知 ， 当 卡 尔 曼 滤 波 增 益 在 k=8 时 已 基 本 达 到 其 稳 定 值018.0019.0k可 写 出 稳 态 时 的 卡 尔 曼 滤 波 器 的 状 态 估 计 方 程 ：)()(833.0090.0361.0627.0)1()1(2121kxkxkxkx)()(833.0091.0361.0627.001)1(248.0200.021kxkxky（ 5-7）或写成)()(833.0090.0361.0627.0) 1() 1(2121kxkxkxkx)(361.0)(627.0)1(248.0200.021kxkxky（5-8）显

8、 然 ， 只 要 选 定 对 系 统 状 态 的 初 始 估 计 值)0(x， 即 可 利 用 依 次 求 得 的 滤波 增 益 值)(kK和 滤 波 估 计 方 程（ 5-3 ）式 对 系 统 状 态 进 行 递 推 估 计 。但 是 ，由于 在 k=10后)(kK值 基 本 趋 于 稳 定 ,故 也 可 利 用 稳 态 滤 波 估 计 方 程（ 5-8 ）式 来对 系 统 状 态 进 行 递 推 估 计 。 此 时 ， 由 于 在 最 初 几 次 测 量 中 采 用 了 较 高 的 滤 波增 益 稳 态 值 ， 对 系 统 状 态 的 前 几 次 估 计 会 出 现 较 大 的 误 差 。

9、 当 然 ， 大 约 在 十个 取 样 周 期 之 后 ， 这 种 情 况 即 可 消 除 。将 上 述 结 果 通 过 MATLAB运 行 之 后 的 图 如 图 5.5-5.8所 示(程 序 见 附 录 ) ：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 12 页012345678910-50510152025Time (sec)Temperature(c)Figure 1 - Temperature (True, Measured, and Estimated)图 5-5 真 实 值 、 测 量 值 及 预 测 估 计 值 随

10、时 间 变 化 的 图在 图 5-5 中 ：红 色 -温 度 的 真 实 值 。蓝 色 -温 度 的 预 测 估 计 值 。绿 色 -温 度 的 测 量 值 。由 图5-5可 看 出 真 实 值 与 预 测 估 计 值 变 化 基 本 一 致 ， 可 看 成 近 似 相 等 ，而 测 量 值 因 受 到 量 测 噪 声 影 响 而 有 所 波 动 变 化 ， 不 过 取 值 在 真 实 值 的 小 范 围内 上 下 波 动 ， 对 以 后 的 预 测 估 计 影 响 较 小 。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 12 页01

11、2345678910-1-0.500.511.5Time (sec)TemperatureError(c)Figure 2 - Temperature Measurement Error and Temperature Estimation Error图 5-6 测 量 值 与 预 测 估 计 值 的 误 差图 5-6将 图 5-5中 的 预 测 估 计 值 与 测 量 值 的 波 动 曲 线 进 行 了 放 大 ， 由 此 ，我 们 可 以 看 出 ， 预 测 估 计 值 的 误 差 基 本 上 趋 向 于 零 ， 说 明 预 测 估 计 值 与 真 实值 趋 于 相 等 ， 可 见 卡 尔

12、 曼 滤 波 器 预 测 估 计 的 准 确 性 。而 测 量 值 误 差 由 于 受 到 量 测 噪 声 的 影 响 则 在 -11摄 氏 度 之 间 上 下 波 动 变化 ， 波 动 幅 度 较 小 ， 对 测 量 结 果 影 响 不 大 。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 12 页012345678910-50510152025Time (sec)Temperature(c/sec)Figure 3 - Temperature (True and Estimated)图 5-7 真实值与预测估计值的变化曲线图 5-7

13、 中说明了真实值与预测估计值的变化趋于一致，近似相等。012345678910-0.15-0.1-0.0500.050.1Time (sec)TemperatureError(c/sec)Figure 4 -Temperature Estimation Error图 5-8 温度预测估计值误差的波动曲线精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 12 页从图5-8 中我们可以看到预测估计值误差的具体波动情况，它的波动变化一直保持在-0.120.1 摄氏度之间，从温度预测估计的角度来说，此误差对整体的预测估计值基本没有影响，从而使得预

14、测估计值与真实值时间的误差基本趋于零。图 5-9 卡 尔 曼 滤 波 增 益 变 化 图由 图 5-9可 知 ， 当 k=8 时 ， 增 益 K(k)逐 渐 趋 向 于 稳 定 。 蓝 线 是 温 度)(1tx的 增益 ， 绿 线 是 温 度)(2tx的 增 益附 录function kalman(duration, dt) % function kalman(duration, dt) % Kalman filter simulation for a vehicle travelling along a road. % INPUTS 精选学习资料 - - - - - - - - - 名师归纳

15、总结 - - - - - - -第 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 ma

16、trix 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 posha

17、t = ; % 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

18、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 Innovatio

19、n 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 * inv(s);% Update the state estimate. xhat = xhat + K * Inn; % Compute the covariance of the estimation error. P = a * P * a - a * P * c * inv(s) * c * P * a + Sw

20、;% Save some parameters for plotting later. pos = pos; x(1); 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 12 页posmeas = posmeas; y; poshat = poshat; xhat(1); vel = vel; x(2); velhat = velhat; xhat(2); end % Plot the results close all; t = 0 : dt : duration; figure; plot(t,pos, t,posmeas, t,

21、poshat); grid; xlabel( Time (sec) );ylabel( Position (feet);title(Figure 1 - Vehicle Position (True, Measured, and Estimated) )figure; plot(t,pos-posmeas, t,pos-poshat); grid; xlabel( Time (sec) );ylabel( Position Error (feet);title(Figure 2 - Position Measurement Error and Position Estimation Error

22、);精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 12 页figure; plot(t,vel, t,velhat); grid; xlabel( Time (sec) );ylabel( Velocity (feet/sec);title(Figure 3 - Velocity (True and Estimated);figure; plot(t,vel-velhat); grid; xlabel( Time (sec) );ylabel( Velocity Error (feet/sec);title(Figure 4 - Velocity Estimation Error);精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 12 页