控制双罐系统
此示例显示如何使用强大的控制工具箱™设计强大的控制器（使用D-K迭代）并对过程控制问题进行稳健性分析。在我们的示例中，该工厂是一个简单的双罐系统。
史密斯等人描述了与该系统有关的其他实验工作。在以下参考文献中：
史密斯，R.S.，J.Doyle，M. Morari和A. Skjellum，“使用MU的案例研究：实验室过程控制问题”，第十二届IFAC世界大会，Vol。8，pp.403-415，1987。
史密斯，R.S和J.Doyle，“两辆坦克实验：基准控制问题”，在“诉讼”美国控制会议上，Vol。3，pp。403-415,1988。
Smith, R.S, and j.c. Doyle，“鲁棒控制模型的不确定性界的闭环继电器估计”，载于第12届IFAC世界大会论文集，第9卷，第57-60页，1993年7月。
植物描述
我们尤其是级联的两个水箱组成，如图1所示。上罐（罐1）通过电机控制的阀通过热和冷水供给。下罐（罐2）通过水从罐1底部的出口供给水。溢流保持罐中的恒定水平。冷水偏置流也馈送罐2并使罐具有不同的稳态温度。
我们的设计目标是控制罐1和2的温度。控制器可以访问参考命令和温度测量。
图1：双罐系统的示意图
坦克变量
让我们给出植物变量以下名称：
FHC.：命令热流执行器FH.:热水流入1号槽FCC.：命令冷流动执行器FC.:冷水流入1号槽F1:从槽1流出的总流量A1：坦克1的横截面积H1.：坦克1水位T1.：坦克的温度1T2.：坦克2的温度A2：坦克2的横截面积H2.：坦克2水位FB.:槽2偏流的流量TB.:槽2偏流温度TH.：热水供应温度TC.：冷水温度
为方便起见，我们定义了一个标准化单位系统，如下所示：
可变单位名称0表示：1表示：-------------------------------温度调节冷水温度。热水温度。高度启动坦克空坦克全流量零零输入流动最大。输入流量
使用上述单位，这些是工厂参数：
A1 = 0.0256;坦克1的％面积（猎杀^ 2）A2 = 0.0477;坦克2的％面积（猎杀^ 2）H2 = 0.241;%水箱2高度，由溢流固定(hunits)fb = 3.28e-5;％偏置流流（猎物^ 3 / sec）FS = 0.00028;％流量缩放（猎杀^ 3 /秒/焦点）th = 1.0;％热水供应温度（调节）Tc = 0.0;%冷水供应温度(单位)结核病= tc;冷偏压流温度(tunits)alpha = 4876;%流量/高度关系常数(hunits/funits)β= 0.59;流量/高度关系（猎物）％常数
变量FS是将输入（0到1个娱乐场）转换为在猎场^ 3 /秒的流量缩放因子。常数alpha和beta描述了罐1的流量/高度关系：
H1 = alpha * f1-beta。
标称坦克模型
我们可以通过在以下操作点周围线性化（所有标准化值）来获得标称罐模型：
H1SS = 0.75;% 1号水箱的水位T1SS = 0.75;% 1号槽温度F1SS =（H1SS + Beta）/ alpha;%流量槽1 ->槽2fs = (th, tc, 1, 1) \ [t1ss * f1ss; f1ss];fhss = fs (1);%热流FCSS = FSS（2）;％冷流T2SS =（F1SS * T1SS + FB * TB）/（F1SS + FB）;% 2号槽温度
一号储罐的标称模型有输入[FH.；FC.[输出[H1.；T1.]：
a = [-1 /（a1 * alpha），0;（beta * t1ss）/（a1 * h1ss）， - （h1ss + beta）/（alpha * a1 * h1ss）];b = fs * [1 /（a1 * alpha），1 /（a1 * alpha）;th / a1，tc / a1];c =α，0;-alpha * T1SS / H1SS，1 / H1SS];d =零（2,2）;tank1nom = ss（a，b，c，d，'InputName', {'fh'那“俱乐部”}，'outputname', {'h1'那't1'}）;CLF步骤（Tank1nom），标题（'坦克1'的阶梯响应）
图2：坦克1的步骤响应。
罐2的标称模型具有输入[| H1 |; | T1 |]和输出T2.：
A = -(h1ss + β + α *fb)/(A2*h2* α);B = [(t2ss + t1ss) /(α* A2 * h2), (h1ss +β)/(α* A2 * h2)];C = 1;D = 0(1、2);tank2nom = ss (A, B, C, D,'InputName', {'h1'那't1'}，'outputname'那't2'）;步骤（TANK2NOM），标题（'坦克2'的步骤响应）
图3:2号罐的阶跃响应。
致动器模型
与执行器相关的动态和饱和性有关，因此我们希望包括执行器模型。在我们使用的频率范围内，我们可以将执行器塑造为具有速率和幅度饱和度的第一订单系统。它是速率限制，而不是极点位置，这限制了大多数信号的执行器性能。对于线性模型，速率限制的一些效果可以包括在扰动模型中。
我们首先使用一个输入（命令信号）和两个输出（驱动信号及其导数）设置执行器模型。我们将在设计控制法时限制致动率时使用衍生品输出。
act_BW = 20;％执行器带宽（RAD / SEC）致动器= [tf(act_BW，[1 act_BW]);tf([act_BW 0]，[1 act_BW])];致动器。OutputName = {'流动'那'流量'};Bodemag（执行器）标题（阀门执行机构动力学的）Hot_act =执行器;设置（hot_act，'InputName'那'fhc'那'outputname', {'fh'那'fh_rate'}）;cold_act =执行器;set（cold_act，'InputName'那'FCC'那'outputname', {“俱乐部”那'fc_rate'}）;
图4：阀门执行机构动力学。
抗锯齿过滤器
所有测量的信号都用四阶巴特沃斯滤波器过滤，每个截止频率为2.25Hz。
FBW = 2.25;%抗混叠滤波器截止(Hz)过滤器= MKFilter（FBW，4，'butterw'）;H1F =过滤器;t1f =过滤器;T2F =过滤器;
模型动态的不确定性
开环实验揭示了系统响应中的一些可变性，并表明线性模型良好良好。如果我们在设计期间未能考虑到此信息，我们的控制器可能会在真实系统上表现不佳。出于这个原因，我们将建立一个不确定性模型，尽可能地密切地匹配我们对物理系统的不确定性的估计。因为模型不确定性或变异的量通常取决于频率，所以我们的不确定性模型涉及频率相关的加权函数来跨频率归一化建模错误。
例如，开环实验表明了大量的动态不确定性T1.回复。这主要是由于混合和热量损失。我们可以将其模拟为乘法（相对）模型错误delta2T1.输出。同样，我们可以将乘法模型错误delta1和delta3添加到H1.和T2.输出如图5所示。
图5:扰动的线性双罐系统的示意图。
为了完成不确定性模型，我们量化了建模错误的函数有多大。虽然难以确定系统中的不确定性的数量，但我们可以根据线性模型准确或差的频率范围寻找粗糙的界限，如在这些情况下：
名义模型
H1.非常准确到至少0.3 Hz。极限环实验
T1.循环表明，不确定性应支配在0.02 Hz以上。有大约180度的额外相位滞后
T1.型号为0.02 Hz。此频率也存在显着的增益损失。这些效果由未确定的混合动态产生。极限环实验
T2.循环表明不确定性应支配在0.03 Hz以上。
此数据表明跨依赖性建模错误界限的以下选择。
WH1 = 0.01 + TF（[0.5,0]，[0.25,1]）;WT1 = 0.1 + TF（[20 * H1S，0]，[0.2,1]）;WT2 = 0.1 + TF（[100,0]，[1,21]）;CLF Bodemag（WH1，WT1，WT2），标题（'建模错误的相对范围'） 传奇（'h1动态'那't1动态'那't2动态'那'地点'那'西北'）
图6:建模错误的相对范围。
现在，我们准备建立不确定的坦克模型，以捕获上面讨论的建模错误。
％标准化错误动态delta1 = ultidyn ('delta1'[1]);delta2 = ultidyn ('delta2'[1]);delta3 = ultidyn（“delta3”[1]);% h1, t1, t2动力学中频率依赖的变异性varh1 = 1 + delta1 * wh1;VART1 = 1 + DELTA2 * WT1;VART2 = 1 + DELTA3 * WT2;％为标称车型添加可变性tank1u =附加（varh1，vart1）* tank1nom;tank2u = Vart2 * Tank2nom;tank1and2u = [0 1;tank2u] * tank1u;
接下来，我们随机采样不确定性，看看建模错误如何影响坦克响应
步骤（TANK1U，1000），标题（'由于建模误差（坦克1）'响应的变异性）
图7:由于建模错误导致的响应变化(坦克1)。
设置控制器设计
现在让我们来看看控制设计问题。我们有兴趣跟踪SetPoint命令T1.和T2.．为了利用h -∞设计算法，我们必须将设计写成闭环增益最小化问题。为此，我们选择捕获扰动特性和性能要求的权重函数，以帮助对相应的频率相关增益约束进行归一化。
这是两个罐问题的合适加权开环传递函数：
图8:双罐系统控制设计互连。
接下来，我们为传感器噪声，设定点命令，跟踪误差和热/冷执行器选择重量。
相对于系统其余部分的动态，传感器动态是微不足道的。传感器噪声不是真的。潜在的噪声源包括热电偶补偿器，放大器和滤波器中的电子噪音，辐射来自搅拌器的辐射噪声，并且接地差。我们使用平滑的FFT分析来估计噪声水平，这表明以下重量：
Wh1noise = zpk (0.01);% h1噪声权重wt1noise = zpk（0.03）;％T1噪音重量wt2noise = zpk（0.03）;% t2噪声权重
错误权重惩罚设定点跟踪错误T1.和T2.．我们将为这些重量选择一流的低通滤波器。我们使用更高的重量（更好的跟踪）T1.因为物理考虑导致我们相信T1.更容易控制而不是T2.．
wt1perf = tf（100，[400,1]）;％t1跟踪误差重量wt2perf = tf（50，[800,1]）;% t2跟踪误差权重CLF Bodemag（WT1Perf，WT2Perf）标题（“频率相关的惩罚在设定点跟踪错误”） 传奇（'t1'那't2'）
图9:对设定值跟踪误差的频率相关惩罚。
参考（设定值）权重反映了这种命令的频率内容。因为流入罐2的大部分水都来自坦克1，改变T2.以变化为主T1.．也T2.通常命令的值接近于T1.．因此，在使用方面表达设定值加权，它更有意义T1.和T2-T1：
t1cmd = wt1cmd * w1
t2cmd = wt1cmd * w1 + wtdiffcmd * w2
其中W1，W2是白噪声输入。足够的体重选择是：
wt1cmd = zpk（0.1）;% t1输入命令权重wtdiffcmd = zpk（0.01）;％T2 - T1输入命令重量
最后，我们希望惩罚执行器的幅度和速率。我们通过加权来这样做FHC.（和FCC.）具有在高频上卷起的功能。或者，我们可以创建一个执行器模型FH.和D | FH | / DT作为输出，并以恒定的重量分开输出。该方法具有减少加权开环模型中的状态数量的优点。
粉虱= ZPK（0.01）;％热执行器罚款WCACT = ZPK（0.01）;％冷执行器罚款Whrate = zpk (50);％热执行器率罚款wcrate = zpk（50）;％冷执行器率罚款
构建加权开环模型
现在我们已经建模了所有工厂组件并选择了我们的设计权重，我们将使用连接函数构建图8所示加权开环模型的不确定模型。
输入= {'t1cmd'那'tdiffcmd'那't1noise'那't2noise'那'fhc'那'FCC'};输出= {'y_wt1perf'那'y_wt2perf'那'是的'那“y_Wcact”那......'y_whrate'那'y_wrate'那'y_wt1cmd'那'y_t1diffcmd'那......'y_t1fn'那'y_t2fn'};hot_act.inputname =.'fhc'；hot_act。OutputName = {'fh''fh_rate'};cold_act。InputName ='FCC'；cold_act。OutputName = {“俱乐部”'fc_rate'};tank1and2u.inputname = {'fh'那“俱乐部”};tank1and2u。OutputName = {'t1'那't2'};t1f.inputname =.'t1'；t1f.outputname ='y_t1f'；t2f.inputname =.'t2'；t2f.outputname ='y_t2f'；Wt1cmd。InputName ='t1cmd'；Wt1cmd。OutputName ='y_wt1cmd'；wtdiffcmd.inputname =.'tdiffcmd'；Wtdiffcmd。OutputName ='y_wtdiffcmd'；粉虱.Inputname =.'fh'；粉虱.outputname =.'是的'；wcact.InputName =.“俱乐部”；wcact.outputname =.“y_Wcact”；whrate.inputname =.'fh_rate'；whrate.outputname =.'y_whrate'；wcrate.Inputname =.'fc_rate'；Wcrate。OutputName ='y_wrate'；wt1perf.inputname =.'u_wt1perf'；wt1perf.outputname =.'y_wt1perf'；Wt2perf。InputName =“u_Wt2perf”；Wt2perf。OutputName ='y_wt2perf'；wt1noise.inputname =.'t1noise'；wt1noise.outputname =.“y_Wt1noise”；wt2noise.inputname =.'t2noise'；wt2noise.outputname =.“y_Wt2noise”；sum1 = sumblk ('y_t1diffcmd = y_wt1cmd + y_wtdiffcmd'）;sum2 = sumblk（'y_t1fn = y_t1f + y_wt1noise'）;sum3 = sumblk（'y_t2fn = y_t2f + y_wt2noise'）;sum4 = sumblk ('u_wt1perf = y_wt1cmd - t1'）;sum5 = sumblk（'u_Wt2perf = y_Wtdiffcmd + y_Wt1cmd - t2'）;％这产生了不确定的状态空间模型P =连接(tank1and2u hot_act、cold_act t1F, t2F, Wt1cmd, Wtdiffcmd, Whact,......Wcact、Whrate Wcrate, Wt1perf、Wt2perf Wt1noise, Wt2noise,......sum1、sum2 sum3、sum4 sum5,输入、输出);disp (加权开环模型:）p
加权开环模型：P =具有10个输出的不确定连续时间 - 空间模型，6个输入，18个状态。模型不确定性包括以下块：Δ1：不确定1x1 lti，峰值增益= 1，1出现Δ2：不确定的1x1 lti，峰值增益= 1，1出现Δ3：不确定1x1 lti，峰值增益= 1，1次出现类型“p.ominalValue“看到标称值”获得（p）“以查看所有属性，以及”p.uncertainty“与不确定元素互动。
H-Infinity控制器设计
通过构造图8的权重和加权开环，我们将控制问题重新循环增益最小化。现在我们可以轻松计算标称坦克模型的增益最小化控制法：
nmeas = 4;%测量次数nctrls = 2;％控件数量[k0，g0，gamma0] = hinfsyn（p.nominalvalue，nmeas，nctrls）;伽玛0.
gamma0 = 0.9016.
可实现的最小闭环增益约为0.9，这向我们表示我们的频域跟踪性能规范由控制器满足k0．在时域中模拟这种设计是检查我们是否正确设置性能权重的合理方法。首先，我们创建一个映射输入信号的闭环模型[T1ref.；T2REF.；t1noise；t2noise.]到输出信号[H1.；T1.；T2.；FHC.；FCC.]：
输入= {'t1ref'那't2ref'那't1noise'那't2noise'那'fhc'那'FCC'};输出= {“y_tank1”那'y_tank2'那'fhc'那'FCC'那'y_t1ref'那'y_t2ref'那......'y_t1fn'那'y_t2fn'};hot_act（1）.inputname ='fhc'；hot_act（1）.outputname ='y_hot_act'；cold_act（1）.inputname ='FCC'；cold_act（1）.outputname ='y_cold_act'；tank1nom。InputName = [hot_act(1)。OutputName cold_act (1) .OutputName];tank1nom。OutputName =“y_tank1”；tank2nom.InputName = tank1nom.outputname;tank2nom.outputname =.'y_tank2'；t1f.inputname = tank1nom.outputname（2）;t1f.outputname ='y_t1f'；t2f.inputname = tank2nom.outputname;t2f.outputname ='y_t2f'；i_tref = zpk（眼睛（2））;i_tref.InputName = {'t1ref'那't2ref'};i_tref.outputname = {'y_t1ref'那'y_t2ref'};sum1 = sumblk ('y_t1fn = y_t1f + t1noise'）;sum2 = sumblk（'y_t2fn = y_t2f + t2noise'）;SIMLFT = CONNECT（TANK1NOM，TANK2NOM，HOT_ACT（1），COLD_ACT（1），T1F，T2F，I_TREF，SUM1，SUM2，输入，输出）;％关闭H-Infinity Controller | K0 | K0 |sim_k0 = lft（simlft，k0）;sim_k0.inputname = {'t1ref'；'t2ref'；'t1noise'；'t2noise'};sim_k0。OutputName = {'h1'；'t1'；'t2'；'fhc'；'FCC'};
现在我们模拟下的闭环响应时的设定值T1.和T2.在80秒和100秒之间：
时间= 0：800;t1ref =（时间> = 80＆time <100）。*（time-80）* - 0.18 / 20 +......（时间> = 100）* - 0.18;t2ref =（时间> = 80和时间<100）。*（时间-80）* - 0.2 / 20 +......（时间> = 100）* - 0.2;t1noise = wt1noise.k * randn（大小（时间））;t2noise = wt2noise.k * randn（大小（时间））;y = lsim（sim_k0，[t1ref; t2ref; t1noise; t2noise]，时间）;
接下来，我们将模拟输出添加到它们的稳态值中，并绘制响应图:
H1 = H1SS + Y（：，1）;t1 = t1ss + y（:,2）;t2 = t2ss + y（:,3）;FHC = FHSS / FS + Y（:,4）;％音符缩放到执行器FCC = FCSS / FS + Y（:,5）;％限制（0 <= FHC <= 1）等
在此代码中，我们绘制了输出，T1.和T2.以及高度H1.坦克1：
情节（时间，H1，“——”、时间t1,' - '、时间、t2,“-”。）;Xlabel（'时间（秒）') ylabel ('测量'） 标题（'H-Infinity Controller K0'的步骤响应） 传奇（'h1'那't1'那't2'）;网格
图10：h -∞控制器k0的阶跃响应。
接下来，我们将命令绘制到冷热执行器。
情节（时间，FHC，' - '、时间、fcc、“-”。）;Xlabel（'时间：秒') ylabel ('执行器'） 标题（h -∞控制器k0的执行器命令） 传奇（'fhc'那'FCC'）;网格
图11：H-Infinity Controller K0的执行器命令。
h -∞控制器的鲁棒性
h∞控制器k0专为标称罐式模型而设计。让我们来看看模型不确定性范围内的扰动模型的票价的票价。我们可以比较标称闭环性能伽玛0.在模型不确定性集上的最坏情况。(更多信息请参见“模型动力学的不确定性”。)
CLPK0 = LFT（P，K0）;％compute和plot最坏情况增益wcsigmaplot (clpk0{1的军医,1 e2}) ylim (-20 [10])
图12：控制器k0的性能分析。
闭环的最坏情况性能比标称性能更差，告诉我们H-Infinity控制器k0对建模错误并不强大。
MU控制器合成
为了解决这种缺乏稳健性，我们将使用斯文设计一个控制器，考虑建模的不确定性，并提供一致的性能为标称和摄动模型。
[kmu，bnd] = musyn（p，nmeas，nctrls）;
D-K ITERATION SUMMARY: ----------------------------------------------------------------- Robust performance Fit order ----------------------------------------------------------------- Iter K Step Peak MU D Fit D 1 8.814 2.862 2.888 4 2 2.504 2.1 2.12 10 3 1.353 1.331 1.339 10 4 1.2 1.198 1.209 18 5 1.194 1.191 1.199 18 6 1.192 1.19 1.195 20 Best achieved robust performance: 1.19
如前所述，我们可以使用控制器模拟闭环响应kmu.
sim_kmu = lft（simlft，kmu）;y = lsim（sim_kmu，[t1ref; t2ref; t1noise; t2noise]，时间）;H1 = H1SS + Y（：，1）;t1 = t1ss + y（:,2）;t2 = t2ss + y（:,3）;FHC = FHSS / FS + Y（:,4）;％音符缩放到执行器FCC = FCSS / FS + Y（:,5）;％限制（0 <= FHC <= 1）等％图| T1 |和| T2 |以及高度| H1 |坦克1.情节（时间，H1，“——”、时间t1,' - '、时间、t2,“-”。）;Xlabel（'时间：秒') ylabel ('测量'） 标题（“MU控制器emu”的阶段响应） 传奇（'h1'那't1'那't2'）;网格
图13：MU控制器KMU的步骤响应。
这些时间响应可与那些k0，并仅显示轻微的性能下降。然而，kmu.更好地对未拼件动态的稳健性更好。
kmu的百分柄性能CLPMU = LFT（P，KMU）;WCSIGMAPLOT（CLPMU，{1E-4,1E2}）ylim（[ - 20 10]）
图14：控制器kmu的性能分析。
您可以使用wcgain直接计算跨频率的最坏情况增益(最坏情况峰值增益或最坏情况h -∞范数)。您还可以计算它对每个不确定元素的灵敏度。结果表明，最坏情况下的峰值增益对的范围变化最为敏感Delta2.．
opt = wcoptions（“敏感”那'在'）;[WCG，WCU，WCINFO] = WCGAIN（CLPMU，OPT）;WCG.
WCG =带有字段的结构：下行：1.3174 Upperbound：1.3201关键频繁：0
wcinfo.sentivity.
ANS =带字段的结构：Delta1：0 Delta2：60 Delta3：10
也可以看看
相关话题
你也可以从以下列表中选择一个网站:
