单容水箱液位控制系统建模与Simulink仿真实践

从一个带本科课设的工程角度讲单容水箱液位控制系统一直是过程控制里最经典的入门被控对象。原因很直接物理过程简单、模型阶次低但又完整覆盖了“对象特性分析—数学模型推导—仿真模型搭建—控制器整定—闭环验证”这一整套工作流。很多同学一开始觉得这块没什么东西一个水箱而已但真正上手做从物料平衡方程到拉氏变换再到Simulink里面跑出闭环响应每一步都有值得细抠的地方。这篇内容我就基于实际做过的仿真工程把单容水箱液位控制系统的模型推导和Simulink仿真完整过一遍写给正在做课设、准备考研复试或者刚进实验室需要搭被控对象模型的读者参考。先说清楚这篇要解决什么问题拿到一个单容水箱对象怎么从物理方程推出传递函数怎么判断它是自衡对象还是非自衡对象又怎么在Simulink里把模型搭出来做PID控制仿真以及搭模型过程中那些容易翻车的细节。内容不追求大而全偏实践和可复现所有参数我都会给出具体的示例值方便对照练习。1. 单容水箱为什么值得单独建模先从被控对象说起1.1 什么是单容对象什么是单容水箱单容对象指的是一类只有一个储能环节的被控对象。对液位系统来说储能环节就是水箱本身它储蓄的是水的体积表现出来的变量是液位高度。所谓单容就是只有一个独立的状态变量描述这个对象只需要一个一阶微分方程。与之相对的双容水箱会有两个串联的水箱两个液位互相影响模型就成了二阶对象后面还会引入纯滞后、耦合等问题复杂度明显上升。单容水箱液位控制系统通常由这几个部分组成进水管路和进水阀或被控泵、水箱本体、出水管路和出水阀或者固定流量的出水泵、液位变送器以及控制器。控制任务一般是维持水箱液位稳定在某个设定值附近或者让液位按一定规律变化。工程上这种场景非常常见比如配料罐液位控制、锅炉汽包水位控制、污水处理中间水池的控制本质上都能抽象成单容或多容液位对象。在这个仿真工程里我的对象定义比较典型一个圆柱形常压水箱底部或侧面有出水口出水口安装了一个手动可调开度的阀门进水侧由一台变频泵或电动调节阀提供可调流量。液位变送器输出标准信号控制器根据液位偏差调节进水流量。整个过程在物理上是一个典型的“进水量-出水量-液位变化”的平衡关系。1.2 模型中的关键假设与边界条件做模型推导之前必须先把假设列清楚这是工程模型和纯数学模型的区别所在。实际水箱的物理过程远比一阶微分方程复杂水面会有波动进水会形成旋涡出水流量和液位之间也不是严格的线性关系。但在课设和初步工程分析阶段我们通常做以下几项假设水箱截面积恒定液位上升或下降过程中有效截面积不变化。对于直筒形的圆柱或方形容器这个假设是成立的但对于顶部缩口或锥形容器就不成立那需要重新推导。液体不可压缩密度恒定不涉及温度变化。水和水的密度差异在这个低压常温场景下可以忽略。进水流量和出水流量不被视为与液位瞬态变化耦合的快速动态。也就是说泵响应速度或阀门动作速度远快于液位变化速度所以在对象模型里不单独考虑执行器动态或者最多用一个一阶惯性环节近似。流出流量与液位高度的关系可以简化成线性关系或者特定的非线性关系。这里需要区分“自衡对象”和“非自衡对象”两种情况它们的建模结果完全不同。自衡对象指的是当进水流量增大导致液位上升时出水流量也会随液位升高而自动增大系统最终会稳定在一个新的液位平衡点上不需要控制器干预也能达到新的稳态。最典型的就是出水口装了一个固定开度的普通阀门出水流量近似与液位高度的平方根成正比。非自衡对象则相反出水流量恒定不随液位变化。比如出水侧用的是定排量泵无论水箱液位多高泵每分钟排出的水量都近似不变。这时候如果进水流量大于出水流量液位就会一直上升直到溢出如果进水流量小于出水流量液位就会一直下降直到抽干。这类对象没有自平衡能力必须依赖控制系统才能稳定运行。这两个边界条件决定了后面的建模方向第一步必须判断清楚做的是哪一种。实际做仿真工程的时候我建议先把这两种对象的模型都写一遍因为很多教材和题库里单容水箱题目会混着出一不小心就把自衡对象的惯性环节当成积分环节去处理结果整个仿真全错了。2. 液位对象传递函数的完整推导从微分方程到拉氏变换2.1 物料平衡方程的建立模型推导的起点是物料平衡方程也就是质量守恒定律在水箱对象上的表达。单位时间内水箱中液体体积的变化量等于流入体积流量减去流出体积流量。用数学公式表达就是A * dh(t)/dt Qin(t) - Qout(t)其中A是水箱截面积单位是平方米h是液位高度单位是米Qin是进水体积流量单位是立方米每秒Qout是出水体积流量单位同样是立方米每秒。这个方程是整个后续推导的根基没有任何近似是对的就是对的只要对象是单容的这个方程就成立。后面要做的所有线性化、拉氏变换、传递函数推导全部围绕这个方程展开。如果只写微分方程当然也能做控制分析但要设计控制器、分析稳定性、整定PID参数工程上还是习惯于把系统表示为传递函数形式。拉氏变换的作用就是把时域微分方程转换成复频域里的代数方程这样一来系统特性就变成输入输出之比这一个简洁的表达式。2.2 线性化处理把这个“根号”拿掉出水流量和水位高度的关系是模型推导里最核心也最容易卡住的地方。如果出水口装的是普通阀门流体力学会告诉我们出水流量的理论公式近似为Qout(t) C * sqrt(h(t))C是一个与阀门开度、出水口面积、流量系数有关的综合系数。也就是说出水流量和液位高度之间是平方根关系不是线性关系。这种情况下整个对象模型的微分方程变成了A * dh(t)/dt Qin(t) - C * sqrt(h(t))这是一个一阶非线性微分方程直接做拉氏变换是做不到的因为拉氏变换要求系统是线性的。所以必须做线性化处理。线性化的思路是在稳态工作点附近把非线性函数展开成一阶泰勒近似。假设系统原来的稳态工作点是(h0, Qin0)此时进出流量相等A * dh/dt 0所以Qin0 C * sqrt(h0)。在工作点附近定义增量变量Δh h - h0ΔQin Qin - Qin0把C * sqrt(h) 在h0附近按一阶泰勒展开忽略高阶无穷小C * sqrt(h) ≈ C * sqrt(h0) (C / (2 * sqrt(h0))) * Δh代入物料平衡方程并把稳态部分抵消掉就可以得到增量形式的线性微分方程A * dΔh/dt ΔQin - (C / (2 * sqrt(h0))) * Δh这里令 R 2 * sqrt(h0) / C它的工程意义是出水阻力代表液位变化单位高度时出水流量的变化率。阻力越大同样的液位变化引起的出流量变化越小系统越不容易自衡。于是方程整理为A * dΔh/dt ΔQin - Δh / R再两边整理成标准的一阶惯性环节形式A * R * dΔh/dt Δh R * ΔQin2.3 自衡对象的传递函数惯性环节对上面整理好的微分方程两边做拉氏变换假设初始条件为零也就是在稳态工作点处开始做增量分析。令T A * RK R则得到T * s * H(s) H(s) K * Qin(s)整理得到传递函数H(s) / Qin(s) K / (T * s 1)这就是典型的一阶惯性环节也叫一阶惯性对象。K是对象的静态增益物理含义是稳态时液位变化量和进水流量变化量之比单位是米每立方米每秒即 s/m²T是时间常数单位是秒物理含义是对象对输入变化的响应快慢程度。T越大液位变化越慢对象越“迟钝”。这里要特别强调一下很多人背公式背得很熟知道一阶惯性环节是K/(Ts1)但不确定K和T各等于什么。实际推导之后一目了然K就是出水阻力RT就是水箱面积A乘以出水阻力R。如果水箱截面积大储水容量大时间常数就大如果出水阀开度小、阻力大时间常数也大。这个对应关系一定要记牢后面整定PID时所有参数都是围绕T和K来算的。我实际用的参数组是这样的水箱直径0.5米截面积A π * 0.25² ≈ 0.196平方米稳态液位h0 0.5米出水阀综合系数C 0.02单位是 m^1.5/s。计算得到稳态出水流量Qout0 0.02 * √0.5 ≈ 0.01414 m³/sR 2 * √0.5 / 0.02 ≈ 70.71 s/m²时间常数T 0.196 * 70.71 ≈ 13.86秒。这个对象的开环响应速度比较慢液位变化要十几秒才能完成大部分过渡过程适合做明显的控制响应演示。2.4 非自衡对象的传递函数积分环节再来看非自衡对象。这种场景下出水流量Qout与液位高度无关是一个恒定值。比如出水侧用的是定流量泵或者出水阀后端有稳压装置总之Qout不随液位变化。物料平衡方程变为A * dh/dt Qin(t) - Qout0其中Qout0是常数。同样在工作点附近做增量分析得到A * dΔh/dt ΔQin(t)两边拉氏变换初始条件为零得到H(s) / Qin(s) 1 / (A * s)这是一个积分环节增益是1/A。如果用执行器输入信号作为输入假设执行器增益为Kv那么传递函数就是Kv/(A*s)等效为一个时间常数为A/Kv的积分对象。积分环节的特点是输入只要有差值输出就一直积累变化不会自动稳定所以这类对象必须靠反馈控制来纠正偏差。自衡和非自衡的模型差异在控制设计上影响非常大。自衡对象可以用比例控制虽然可能有稳态误差但至少不会发散非自衡对象如果只用比例控制增益不合适时会变成等幅振荡甚至发散。所以拿到题目先做对象特性判断比急着仿真要重要得多。2.5 用表格对照两种情况为了方便对照我把两种单容液位对象的模型关键点整理成表格仿真搭模型之前先对上号对比项自衡对象非自衡对象出水侧结构普通阀门出流量随液位变化定流量泵或恒出流装置出流量与液位关系Qout ≈ C * sqrt(h) 或 Qout h/RQout为常数线性化后模型T * dh/dt h R * QinA * dh/dt Qin传递函数K/(Ts1)一阶惯性环节1/(As)积分环节稳态特性存在新的平衡点有自衡能力无平衡点液位持续积分变化控制难度相对容易比例控制也有一定余量必须闭环整定不好易发散这个表格基本把两种对象的本质区别抓出来了。实际做Simulink仿真时我建议两种模型都建一下分别跑开环阶跃响应感受一下曲线的区别。这比对着PPT看十遍都管用。3. Simulink工程搭建把数学模型接成仿真回路3.1 方案选择直接用传递函数还是搭模块Simulink里实现单容水箱对象有两条路。第一条路最省事直接在Transfer Fcn模块里填入K/(Ts1)设置好K和时间常数再给一个阶跃输入一个Scope就能看开环响应。这个方案适合验证传递函数推导是否正确也适合快速跑到闭环做PID整定。第二条路是从微分方程直接搭用积分器、增益、加法器搭出A*dh/dt Qin - h/R这个结构。这个方案能直观看到物理过程中每个变量的走向也方便后续加入对象非线性比如平方根出水特性。我的建议是两条腿走路。先用传递函数模块把闭环控制回路跑通验证控制器参数和整定方法再用积分器搭物理结构模型把线性对象和带平方根的非线性对象做对比观察线性化近似在什么范围有效。实际课设展示时这个对比本身就是很好的加分项。还要提醒一点传递函数模块内部默认初始状态是零。如果你的仿真设定液位从0.5米开始就得在积分器的初始条件里设置或者通过把输入信号叠加上稳态值来模拟工作点。很多同学的模型跑出来曲线不符合物理直觉就是因为初始状态没设对仿真一开始液位从零开始往上爬看起来极其奇怪。3.2 模块连接与参数设置下面给出一个可复现的自衡对象线性模型搭建步骤使用Simulink标准库模块输入信号使用Step模块阶跃时间设为0初始值0终值0.001代表进水流量突增0.001 m³/s。进水流量信号进入一个Gain模块增益为K70.71对应传递函数分子。之后接一个一阶惯性环节但这里更直接的办法是用一个Transfer Fcn模块分子为[70.71]分母为[13.86, 1]。液位输出从Transfer Fcn出来之后接Scope显示同时用To Workspace模块导出时间序列方便后续在MATLAB里做数据分析。若要搭建闭环回路把液位信号反馈给PID ControllerPID输出作为执行器指令执行器输出经过饱和限幅后作为进水流量叠加到对象上。这里需要特别说明的是设定值信号的单位。液位设定值的物理单位是米而PID Controller输出是执行器指令需要经过一个阀门或泵的增益才能换算成流量。如果直接拿PID输出当流量K值就必须包含执行器增益。我建议在模型中显式加一个执行器增益模块比如阀门的增益Kv 0.01 m³/s每单位指令。这样传递关系就变成液位偏差 → PID → 执行器增益 → 流量 → 对象 → 液位。模块参数设置示例模块参数设置说明StepStep time 0Initial value 0Final value 0.001进水流量增量Transfer FcnNumerator [70.71]Denominator [13.86, 1]对象模型 K/(Ts1)PID Controller按整定结果设置见后文控制器SaturationUpper limit 0.002Lower limit 0执行器限幅To WorkspaceVariable name h_outSave format Timeseries导出液位数据搭完模型先别急着加控制器跑一次开环阶跃。理论上液位从0.5米稳态开始经过约三倍时间常数后变化量接近0.0707米最终液位约0.5707米。如果跑出来的结果不是这个数说明K或者T的设置有问题先排查模型再谈控制。3.3 加入执行器和测量噪声的仿真模型实际系统里没有理想模块。电动调节阀有动作惯性液位变送器有测量噪声信号传输有延迟。为了让仿真更接近真实工程可以在基本模型上逐步增加这些环节。执行器动态我用一个一阶惯性环节来近似时间常数设为0.5秒表示阀门从收到指令到流量实际变化需要一定响应时间。这个惯性放在PID输出和对象输入之间。加上执行器惯性后系统从一阶对象变成了二阶系统控制难度会上升一点比例增益太大时会出现超调甚至振荡。测量噪声用Random Number模块模拟均值设为0方差根据传感器精度来设置。比如液位变送器精度为满量程的0.5%量程1米则噪声标准差取0.005米左右。噪声加到液位测量信号上然后进入PID的反馈通道。加上噪声之后就能明显看出一个问题纯微分项对测量噪声非常敏感D项太大会把噪声放大成剧烈的执行器抖动。这是实际工程里最常踩的坑。信号延迟用Transport Delay模块模拟延迟时间设为0.5秒到1秒。液位信号从变送器到控制器的传输滞后虽然时间不长但会让控制系统稳定性变差。加了延迟之后再去整定P就能比较容易复现出临界振荡现象用临界比例度法的时候这个延迟模块几乎是必须的。4. PID参数整定与闭环仿真只看响应曲线说话4.1 开环阶跃响应的数据读取模型搭好之后第一步是做开环阶跃响应实验。我给对象输入一个幅度为0.001 m³/s的阶跃流量增量记录液位响应曲线。从曲线上读取稳态变化量Δh_ss 0.0707米得到过程增益K 0.0707 / 0.001 70.71和理论计算值完全一致。再看响应达到稳态的63.2%所经过的时间约为13.86秒得到时间常数T。这个开环实验有两个作用。其一验证了模型参数和推导的一致性避免带着错误模型去做闭环整定。其二为后面的过程反应曲线整定法提供了原始数据。一阶对象没有纯滞后的时候临界比例度法会失效因为无论比例增益多大一阶闭环系统都不会持续振荡。因此我选择在模型中增加0.5秒的传输延迟再去做整定这样更贴近实际也方便使用经典整定公式。4.2 过程反应曲线法整定过程反应曲线法也叫开环阶跃响应法适合对象传递函数能近似成一阶惯性加纯滞后的情况。我们的对象加上传输延迟后近似为G(s) K * e^(-Ls) / (T * s 1)其中K 70.71T 13.86秒L 0.5秒。从阶跃响应曲线上读到滞后时间L和等效时间常数T之后按齐格勒-尼科尔斯开环整定公式计算控制器类型KpTiTdPT / (K * L)——PI0.9 * T / (K * L)L / 0.3—PID1.2 * T / (K * L)2 * L0.5 * L代入数值计算Kp 1.2 * 13.86 / (70.71 * 0.5) ≈ 0.47。这是一个无量纲比例增益对应到Simulink里的PID模块如果系统输入输出都是工程单位这个值可以直接使用。Ti 2 * 0.5 1秒Td 0.5 * 0.5 0.25秒。这个整定结果往往超调偏大工程上一般会在此基础上适当减小比例增益、增大积分时间让系统更稳健。课设仿真里可以保留两组参数做对比一组是齐格勒-尼科尔斯参数一组是手动软化后的参数分别跑闭环响应用超调量和调节时间说明整定方法的适用性。4.3 在Simulink里做PID参数实验对比把PID Controller模块接入闭环回路之后注意几个设置控制器输出限幅要设成和执行器一致避免PID内部积分项无限制累积PID模块的积分器和微分器初始状态也要合理设置最好让控制作用从零偏置开始。我这里设定值从0.5米阶跃到0.55米观察液位跟踪情况仿真时长设为60秒。分别用手动经验法试三组参数第一组只加比例控制Kp 0.2。液位最终稳定在0.535米附近有大约0.015米的稳态误差。原因很直观比例控制只能缩小误差不能消除误差要想消除稳态误差必须引入积分作用。第二组PI控制Kp 0.3Ki 0.05。稳态误差被彻底消除液位最终准确稳定在0.55米。但响应过程出现了一个小超调大约4%调节时间在25秒左右。这是因为积分项在误差累积过程中产生了过冲之后逐渐回调。第三组PID控制Kp 0.4Ki 0.1Kd 0.1。超调进一步减小响应速度更快调节时间缩短到12秒左右。但如果把Kd继续调大会出现执行器频繁抖动尤其在加了测量噪声的模型里高频噪声被微分放大控制曲线变成毛刺状。把这三种情况的响应曲线放在同一张图里对比能非常清楚地向读者展示P、PI、PID各自的作用。这也是我建议大家在博文或课设报告中必放的图。Simulink里用多个Scope合并或者都用To Workspace导出后用MATLAB的plot命令画在一张图效果都比直接截图Scope好看得多。4.4 为什么PI比P更适合液位控制从工程实践角度说液位控制对象带积分或者近似带积分特性并不少见。即便是自衡对象只要时间常数比较大比例控制的稳态误差也不太容易接受。以我这个对象为例纯比例控制Kp 0.2时稳态误差是0.015米相对设定值变化0.05米来说已经达到30%这在工业液位控制中是不可接受的。液位控制的目标本来就是“稳定在设定值”所以工程上绝大多数液位回路都使用PI控制。微分项在液位控制里的地位比较微妙。液位变化本身是大惯性慢过程微分作用对快速扰动有一定抑制效果但对测量噪声特别敏感。实际工业液位控制器里D参数经常设得很小甚至直接设为零就是因为液位变送器在液面波动时输出噪声偏大微分项容易引发调节阀频繁动作。仿真模型里加不加测量噪声对微分项优缺点的体会完全不一样所以建议大家务必在模型里加上噪声模块再做结论。5. 仿真结果分析与常见的翻车现场5.1 典型阶跃响应曲线的判读闭环仿真的结果不能只看最终稳定不稳定还要看动态过程的几个硬指标超调量、上升时间、峰值时间、调节时间和稳态误差。超调量看的是液位会不会冲破上限甚至溢出调节时间看的是系统受扰动后多久能恢复稳态误差看的是最终精度。一个比较理想的液位PID响应曲线应该是这样的液位从初始值快速上升接近设定值时出现一个轻微超调然后回落并稳定在设定值附近整个过程没有剧烈振荡执行器输出没有频繁满幅。如果响应曲线出现等幅振荡大概率是比例增益过大如果响应缓慢且长时间回不到设定值大概率是积分作用太弱如果曲线前期猛冲、后期拖着长长的尾巴可能是积分饱和或者微分参数不合适。我在仿真中特别注意观察执行器输出曲线。很多时候液位曲线看起来还可以但阀门开度输出一直在上下限之间来回撞这种系统在真实设备上是没法用的。执行器频繁大幅动作会快速磨损阀门所以控制器整定不能只看被控量还要看操纵量的品质。这一点课设报告里很少有人写但面试或实际项目答辩时往往会被问到。5.2 单位换算错误是第一个坑单容水箱仿真翻车最集中的原因就是单位没有统一。教材例题里经常把流量写成升每分钟液位写成厘米但传递函数推导时用的又是米和立方米每秒混在一起之后K和T的数值就差出好几个数量级。我见过有同学Simulink跑出来的液位响应曲线是负数排查半天正是因为进水流量用了立方米每小时而对象增益算的是立方米每秒对应的响应差了3600倍。建模型之前先把所有变量和参数的单位写在纸上液位米流量立方米每秒截面积平方米阀门阻力秒每平方米时间常数秒。所有模块参数都按这套单位填。如果输入信号想用立方米每小时那就在进入对象前用Gain模块做一次3600倍的单位换算不要在地层深处偷偷改数字。另外Simulink的Transfer Fcn模块默认信号是连续的数值量本身没有单位校验。模型不会提醒你单位错了你只能靠开环稳态值来反推。我的检查方法是给一个已知的阶跃流量增量看理论计算的稳态液位增量是否和仿真一致。如果差得远先查单位再查参数。5.3 阀门非线性与积分饱和线性模型里出水流量按h/R处理实际阀门是平方根特性。两者的差异在稳态工作点附近很小但工作点变化大或者液位接近零时差异就明显了。比如液位很低时平方根特性使得出流量变化率很大实际对象比线性模型更“敏感”液位很高时出流量变化率变小对象变得比线性模型更“迟钝”。更麻烦的是执行器饱和。PID输出限幅之后控制器算出的理想控制量可能超过执行器能力此时误差还在积分项继续累积等误差反向时积分值才慢慢释放就会出现明显的超调甚至长时间无法收敛。这个现象叫积分饱和解决的办法是使用带抗积分饱和功能的PID结构。Simulink的PID Controller模块里提供了Anti-windup method选项选成back-calculation再设一个恰当的Back-calculation coefficient就能明显改善。我在非线性模型里故意把PI参数调大触发积分饱和液位曲线出现了幅度很大的第一峰然后在设定值附近来回震荡了好几次才稳定。换成抗积分饱和设置之后第一个超调峰值明显降低后续振荡也更快平息。这个对比实验在课设展示里很有说服力。5.4 变步长求解器与信号采样Simulink默认使用变步长求解器比如ode45对简单线性对象没什么问题。但加了传输延迟、继电器、饱和、PID抗积分饱和这些非线性模块之后变步长求解器在状态切换时刻会大幅减小步长仿真速度变慢甚至可能因为事件检测过多而报错。我建议这类控制仿真统一改成固定步长求解器比如ode4步长取0.01秒到0.02秒。这样做的好处是仿真结果确定性好同一组参数每次跑出来完全一样方便对比实验。步长也不是越小越好。步长太小会成倍增加仿真时间意义不大步长太大则会在快速动态上引入数值误差。对于时间常数13.86秒这个对象0.01秒步长已经足够精细。如果模型中加入了0.5秒的纯延迟最好把延迟模块的缓冲大小设成步长的整数倍避免插值误差。另外要养成好习惯仿真结果不要只在Scope里看用To Workspace导出数据之后用MATLAB脚本统一计算超调量、调节时间等指标。手动读Scope曲线误差很大特别是调节时间往往要看到1秒以内才算准眼睛盯屏幕读不出来。写一个简单脚本用find函数找稳态误差进入2%范围的时间点几行代码搞定报告里放的数据也更有说服力。最后再分享一个实际心得搭这类模型时先跑最简线性闭环再一步步加延迟、加执行器饱和、加测量噪声每一步加完都跑一次看看曲线。很多人一上来就搭一个非常复杂的全功能模型结果出了问题都不知道是哪个环节引起的。分层搭建、逐层验证是Simulink仿真工程里最值得养成的习惯。我每次带人做液位控制仿真都强调这个流程按这个顺序来基本一天之内能把模型、整定和结果分析全部做完而且每个结论都有清晰的仿真依据。