双罐系统的控制
这个例子展示了如何使用鲁棒控制工具箱™设计一个鲁棒控制器(使用D-K迭代)，并对过程控制问题进行鲁棒分析。在我们的例子中，工厂是一个简单的双罐系统。
Smith等人在以下文献中描述了与该系统相关的其他实验工作:
Smith, R.S, J. Doyle, M. Morari和A. Skjellum，“一个使用mu的案例研究:实验室过程控制问题”，第10届IFAC世界大会论文集，第8卷，第403-415页，1987年。
Smith, r.s.和J. Doyle，“双坦克实验:一个基准控制问题”，《美国控制会议论文集》，第3卷，第403-415页，1988年。
Smith, r.s.和J. C. Doyle，“鲁棒控制模型的不确定界的闭环中继估计”，《第十二届IFAC世界大会论文集》，第9卷，第57-60页，1993年7月。
植物的描述
在我们的例子中，工厂由两个串联的水箱组成，如图1所示。上面的水箱(水箱1)通过计算机控制的阀门输入热水和冷水。下槽(槽2)由槽1底部的出口进水。2号槽的溢流保持恒定水平。冷水偏置流也为槽2提供水，使槽2具有不同的稳态温度。
我们的设计目标是控制1号罐和2号罐的温度。控制器可以访问参考命令和温度测量。
图1:双罐系统原理图
坦克变量
让我们给植物变量命名如下:
fhc:热流执行器命令跳频:热水流入1号槽fcc:命令到冷流执行器足球俱乐部:冷水流入1号槽f1:槽1的总流出量A1:槽1截面积h1:槽1水位t1:槽1温度t2:槽2温度A2:槽2截面积h2:槽2水位神奇动物:槽2偏置流流量结核病:槽2偏置流温度th:热水供应温度tc:冷水供水温度
为方便起见，我们将标准化单位系统定义如下:
变量单位名称0表示:1表示:-------- --------- -------- --------温度tunit冷水温度热水温度高度hunit槽空槽满流量funit零输入流量最大。输入流
使用上述单位，这些是工厂参数:
A1 = 0.0256;1号槽面积% (hunits^2)A2 = 0.0477;2号槽面积% (hunits^2)H2 = 0.241;槽2高度，溢出固定(hunits)Fb = 3.28e-5;%偏置流流量(hunits^3/sec)Fs = 0.00028;%流量缩放(hunits^3/秒/funit)Th = 1.0;%热水供应温度(tunits)Tc = 0.0;冷水供应温度(tunits)TB = tc;冷偏置流温度(tunits)Alpha = 4876;%恒定的流量/高度关系(hunits/funits)Beta = 0.59;流量/高度关系%常数(hunits)
变量fs是一个流量缩放因子，它将输入(0到1 funits)转换为hunits^3/秒的流量。常数和描述了1号罐的流量/高度关系:
H1 = *f1-。
标称坦克型号
我们可以通过围绕以下工作点(所有归一化值)线性化得到标称坦克模型:
H1ss = 0.75;1号槽水位%T1ss = 0.75;槽1温度%F1ss = (h1ss+beta)/alpha;%流量罐1 ->罐2FSS = [th,tc;1,1]\[t1ss*f1ss;f1ss];FHSS = fss(1);%热流量FCSS = fss(2);%冷流T2ss = (f1ss*t1ss + fb*tb)/(f1ss + fb);2 .槽温度%
罐体1的公称型号有输入[跳频；足球俱乐部和输出[h1；t1]：
A = [-1/(A1*alpha)， 0;(β* t1ss) / (A1 * h1ss)、- (h1ss +β)/(α* A1 * h1ss)];B = fs*[1/(A1*alpha)， 1/(A1*alpha);th / A1, tc / A1];C = [alpha, 0;t1 * t1ss / h1ss, 1 / h1ss];D = 0 (2,2);(A,B,C,D，“InputName”, {“跳频”，“俱乐部”}，“OutputName”, {“标题”，“t1”});CLF step(tank1nom)， title(“1号罐的阶梯反应”）
图2:1号槽的阶跃响应。
罐体2的标称模型有输入[|h1|;|t1|]和输出t2：
A = -(h1ss + beta + alpha*fb)/(A2*h2*alpha);B = [(t2ss + t1ss) /(α* A2 * h2), (h1ss +β)/(α* A2 * h2)];C = 1;D = 0 (1,2);tank2nom = ss(A,B,C,D，“InputName”, {“标题”，“t1”}，“OutputName”，《终结者2》）;步骤(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_rate”});cold_act =执行机构;集(cold_act,“InputName”，“fcc”，“OutputName”, {“俱乐部”，“fc_rate”});
图4:阀门执行机构动力学。
抗锯齿过滤器
所有测量信号都用四阶巴特沃斯滤波器进行滤波，每个滤波器的截止频率为2.25 Hz。
FBW = 2.25;%抗混叠滤波器截止(Hz)Filter = mkfilter(fbw,4，“Butterw”）;h1F =过滤器;t1F = filter;t2F = filter;
模型动力学的不确定性
开环实验表明系统响应具有一定的可变性，线性模型在低频时表现良好。如果我们在设计过程中没有考虑到这些信息，那么我们的控制器在实际系统中可能表现不佳。因此，我们将建立一个不确定性模型，尽可能接近我们对物理系统不确定性的估计。由于模型不确定性或可变性的数量通常取决于频率，我们的不确定性模型涉及频率相关的加权函数，以跨频率归一化建模误差。
例如，开环实验表明在t1响应。这主要是由于混合和热损失。我们可以将其建模为乘(相对)模型误差在t1输出。类似地，我们可以将模型误差Delta1和Delta3添加到h1而且t2输出如图5所示。
图5:摄动线性双罐系统的示意图。
为了完成不确定性模型，我们量化了建模误差的大小作为频率的函数。虽然很难精确地确定系统中的不确定性的数量，但我们可以根据线性模型准确或糟糕的频率范围寻找粗略的边界，如以下情况:
的标称模型
h1是非常精确的，至少可以达到0.3 Hz。极限环实验
t1循环显示不确定度应在0.02 Hz以上占主导地位。有大约180度的附加相位滞后
t1模型的频率约0.02赫兹。在这个频率上也有显著的增益损失。这些效应来自于未建模的混合动力学。极限环实验
t2循环显示不确定度应在0.03 Hz以上占主导地位。
这些数据为频率相关的建模误差边界提供了以下选择。
Wh1 = 0.01+tf([0.5,0]，[0.25,1]);Wt1 = 0.1+tf([20*h1ss,0]，[0.2,1]);Wt2 = 0.1+tf([100,0]，[1,21]);clf bodemag(Wh1,Wt1,Wt2)， title(“建模误差的相对边界”)传说(“h1动力学”，“t1动力学”，“t2动力学”，“位置”，“西北”）
图6:建模误差的相对边界。
现在，我们准备构建不确定的坦克模型，捕捉上面讨论的建模错误。
规范化误差动态Delta1 = ultidyn(“delta1”[1]);delta = ultidyn(“delta2”[1]);delta = ultidyn(“delta3”[1]);% h1, t1, t2动态中频率相关的变异性varh1 = 1+delta1*Wh1;vart1 = 1+delta2*Wt1;vart2 = 1+delta3*Wt2;增加标称模型的可变性Tank1u = append(varh1,vart1)*tank1nom;Tank2u = vart2*tank2nom;Tank1and2u = [0 1;tank2u] * tank1u;
接下来，我们对不确定性进行随机抽样，看看建模误差如何影响坦克的响应
步骤(tank1u, 1000)、标题(“由于建模错误导致的响应的可变性(坦克1号)”）
图7:由于建模错误导致的响应的可变性(Tank 1)。
设置控制器设计
现在让我们着眼于控制设计问题。我们感兴趣的是跟踪设定值命令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) title(“定点跟踪误差的频率依赖性惩罚”)传说(“t1”，《终结者2》）
图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)，带有一个在高频下滚动的函数。或者，我们可以创建一个驱动器模型跳频d|fh|/dt作为输出，分别用恒权对每个输出进行加权。这种方法的优点是减少了加权开环模型的状态数。
Whact = 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_Whact”，“y_Wcact”，.．.“y_Whrate”，“y_Wcrate”，“y_Wt1cmd”，“y_t1diffcmd”，.．.“y_t1Fn”，“y_t2Fn”};hot_act。InputName =“fhc”；hot_act。OutputName = {“跳频”“fh_rate”};cold_act。InputName =“fcc”；cold_act。OutputName = {“俱乐部”“fc_rate”};tank1and2u。InputName ={“跳频”，“俱乐部”};tank1and2u。OutputName = {“t1”，《终结者2》};t1F。InputName =“t1”；t1F。OutputName =“y_t1F”；t2F。InputName =《终结者2》；t2F。OutputName =“y_t2F”；Wt1cmd。InputName =“t1cmd”；Wt1cmd。OutputName =“y_Wt1cmd”；Wtdiffcmd。InputName =“tdiffcmd”；Wtdiffcmd。OutputName =“y_Wtdiffcmd”；Whact。InputName =“跳频”；Whact。OutputName =“y_Whact”；Wcact。InputName =“俱乐部”；Wcact。OutputName =“y_Wcact”；Whrate。InputName =“fh_rate”；Whrate。OutputName =“y_Whrate”；Wcrate。InputName =“fc_rate”；Wcrate。OutputName =“y_Wcrate”；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 = connect(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状态。模型不确定度由以下块组成:delta1:不确定1x1 LTI，峰值增益= 1,1次delta2:不确定1x1 LTI，峰值增益= 1,1次delta3:不确定1x1 LTI，峰值增益= 1,1次键入"P. nominalvalue "查看标称值，"get(P)"查看所有属性，以及"P. uncertainty "与不确定元素交互。
h -∞控制器设计
通过构造图8的权重和权重开环，我们将控制问题重新定义为闭环增益最小化问题。现在我们可以很容易地计算标称坦克模型的增益最小控制律:
Nmeas = 4;测量次数%NCTRLS = 2;%控件数量[k0,g0,gamma0] = hinsyn (P.NominalValue,nmeas,nctrls);gamma0
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). 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(eye(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,input,outputs);以h -∞控制器|k0|闭合环路Sim_k0 = lft(simlft,k0);sim_k0。InputName = {“t1ref”；“t2ref”；“t1noise”；“t2noise”};sim_k0。OutputName = {“标题”；“t1”；《终结者2》；“fhc”；“fcc”};
现在我们模拟了在降低设定值时的闭环响应t1而且t280秒到100秒之间:
时间= 0:800;t1ref =(> = 80 &时间< 100)。*(时间- 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以及身高h11号槽:
情节(h1,“——”、时间t1,“- - -”、时间、t2,“-”。）;包含(的时间(秒)) ylabel (“测量”)标题(h -∞控制器k0的阶跃响应)传说(“标题”，“t1”，《终结者2》）;网格
图10:h∞控制器k0的阶跃响应。
接下来，我们绘制到热和冷执行器的命令。
情节(时间、fhc“- - -”、时间、fcc、“-”。）;包含(“时间:秒”) ylabel (“执行机构”)标题(h -∞控制器k0的执行器命令)传说(“fhc”，“fcc”）;网格
图11:h -∞控制器k0的执行器命令。
h -∞控制器的鲁棒性
h∞控制器k0是为标称型号的坦克而设计的。让我们看看它在模型不确定性边界内的摄动模型的效果如何。我们可以比较标称闭环性能gamma0最坏情况下的性能超过模型的不确定性集。(更多信息请参见“模型动力学的不确定性”。)
clpk0 = lft(P,k0);计算并绘制最坏情况增益图。wcsigmapplot (clpk0，{1e-4,1e2}) ylim([-20 10])
图12:控制器k0的性能分析。
闭环的最坏情况性能明显比标称性能差这告诉我们h∞控制器k0对建模误差不健壮。
Mu控制器综合
为了弥补健壮性的不足，我们将使用musyn设计一种考虑建模不确定性并对名义模型和摄动模型提供一致性能的控制器。
[kmu,bnd] = musyn(P,nmeas,nctrls);
D-K迭代总结 : ----------------------------------------------------------------- 强劲的性能符合订单 ----------------------------------------------------------------- Iter K步峰μD适合D 1 8.814 2.862 2.888 4 2 4 10 3 1.351 1.329 1.337 2.497 2.091 2.112 1.201 1.198 1.209 18 5 18 6 1.194 1.191 1.199 1.192 1.19 1.195 20最好实现健壮的性能:1.19
和以前一样，我们可以用控制器模拟闭环响应kmu
Sim_kmu = lft(simlft,kmu);Y = lsim(sim_kmu，[t1ref; t1ref; t1noise;t2noise]，time);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|以及水箱1的高度|h1|情节(h1,“——”、时间t1,“- - -”、时间、t2,“-”。）;包含(“时间:秒”) ylabel (“测量”)标题(“控制器kmu的阶跃响应”)传说(“标题”，“t1”，《终结者2》）;网格
图13:控制器kmu的阶跃响应。
这些时间反应与那些k0，只显示出轻微的性能下降。然而,kmu对于未建模的动态的健壮性更好。
kmu的最坏情况表现clpmu = lft(P,kmu);wcsigmapplot (clpmu，{1e-4,1e2}) ylim([-20 10])
图14:kmu控制器的性能分析。
你可以用wcgain直接计算跨频率的最差增益(最差峰值增益或最差h∞范数)。你也可以计算它对每个不确定元素的灵敏度。结果表明，最坏情况峰值增益对范围的变化最为敏感delta2．
opt = wcOptions(“敏感”，“上”）;[wcg,wcu,wcinfo] = wcgain(clpmu,opt);wcg
wcg = struct with fields: LowerBound: 1.3174 UpperBound: 1.3201 CriticalFrequency: 1.9350e-11
wcinfo。灵敏度
Ans = struct with fields: delta1: 0 delta2: 60 delta3: 10
另请参阅
相关的话题
您也可以从以下列表中选择网站: