残差分析を使用した遠心ポンプの故障診断

この例では次を使用します:

この例では,ポンプシステムで発生するさまざまなタイプの故障を検出して診断するモデルパリティ方程式に基づく方法を説明します。ここでは定常状態実験を使用した遠心ポンプの故障診断で提示されている手法を，ポンプの動作が複数の動作状態にわたる状況に拡張します。

ここではRolf Isermann著の書籍“故障诊断应用”で紹介されている遠心ポンプ解析を使用します[1]。系统识别工具箱™、Statistics and Machine Learning Toolbox™、Control System Toolbox™、および Simulink™ の機能を使用しますが、Predictive Maintenance Toolbox™ は必要ありません。

複数速度のポンプ運転-残差分析による診断

定常状態の揚程とトルクの方程式では,速度が急激に変化したり速度範囲がさらに広くなる状態でポンプを運転した場合,正確な結果を得ることができません。摩擦やその他の損失の影響が増し，モデルのパラメ，タ，が速度に依存するようになる可能性があります。そのような場合に広く応用できる方法として,動作のブラックボックスモデルを作成することができます。そうしたモデルのパラメ，タ，は，物理的に有意である必要はありません。このモデルは既知の動作のシミュレションのためのデバスとして使用されます。モデルの出力を対応する測定信号から差し引いて"残差"を計算します。平均，分散およびべき乗などの残差のプロパティは，通常動作と故障動作を区別するために使用されます。

静的な揚程の方程式および動的なポンプ——配管の方程式を使用して,図に示されている4つの残差を計算することができます。

このモデルには次のコンポ，ネントがあります。

静的なポンプモデル: $Δ p_{}^{ˆ} （ t ）＝ θ_{1} ω^{2} （ t ） + θ_{2} ω （ t ）$
動的な配管モデル: $问_{}^{ˆ} （ t ）＝ θ_{3.} + θ_{4} \sqrt{Δ p （ t ）} + θ_{5} 问_{}^{ˆ} （ t - 1 ）$
動的なポンプ-配管モデル: ${问_{}^{ˆ}}_{}^{ˆ} （ t ）＝ θ_{3.} + θ_{4} \sqrt{{Δ p （ t ）}_{}^{ˆ}} + θ_{5} {问_{}^{ˆ}}_{}^{ˆ} （ t - 1 ）$
動的な逆ポンプモデル: ${米_{}^{ˆ}}_{米 o t o r} （ t ）＝ θ_{6} + θ_{7} ω （ t ） + θ_{8} ω （ t - 1 ） + θ_{9} ω^{2} （ t ） + θ_{10} {米_{}^{ˆ}}_{米 o t o r} （ t - 1 ）$

モデルパラメタ $θ_{1} ，．．．， θ_{10}$ はポンプ回転数への依存を示しています。この例では，パラメ，タ，の区分的線形近似を計算します。動作領域を3の領域に分割します。

$ω \leq 900 R P 米$
$900 < ω \leq 1500 R P 米$
$ω > 1500 R P 米$

健全なポンプが,基準速度範囲0 ~ 3000 RPMにおいて,閉ループのコントローラーをもつ閉ループで実行されています。基準入力は修正PRBS信号です。モータートルク,ポンプトルク,速度および圧力の測定値は10 Hzのサンプリング周波数で収集しました。測定した信号を読み込み,基準と実際のポンプ速度をプロットします(これらは最大20 MBの大規模なデータセットで,MathWorksサポートファイルサイトで入手します)。

url =“//www.tatmou.com/金宝appsupportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/DynamicOperationData.mat”；websave (“DynamicOperationData.mat”url);负载DynamicOperationData图(t, RefSpeed, t, w)“时间(s)”) ylabel (“泵转速(RPM)”)传说(“参考”，“实际”）

ポンプ速度の範囲に基づいて動作範囲を定義します。

I1 = w<=900;%第一工作状态I2 = w>900 & w<=1500;%第二工作状态I3 = w>1500;%第三工作状态

モデルの同定

答:静的なポンプモデルの同定

ポンプ回転数 $ω （ t ）$ と圧力差 $Δ p （ t ）$ の測定値を入出力デ，タとして使用して，静的なポンプ方程式のパラメ，タ， $θ_{1}$ および $θ_{2}$ を推定します。この推定を実行する補助関数staticPumpEstを参照してください。

Th1 = 0 (3,1);Th2 = 0 (3,1);dest = nan(大小(dp));估计压差%[th1(1)， th2(1)， dest (I1)] = staticPumpEst(w, dp, I1);区域1的估计值[th1(2)， th2(2)， dest (I2)] = staticPumpEst(w, dp, I2);区域2的估计值[th1(3)， th2(3)， dest (I3)] = staticPumpEst(w, dp, I3);区域3的估计值Plot (t, dp, t, dest)%比较实测和预测的压差包含(“时间(s)”) ylabel (‘\δP’)传说(“测量”，“估计”，“位置”，“最佳”)标题(“静态泵模型验证”）

B.動的な配管モデルの同定

流量 $问（ t ）$ と圧力差 $Δ p （ t ）$ の測定値を入出力デ，タとして使用して，配管吐出流の方程式 $问_{}^{ˆ} （ t ）＝ θ_{3.} + θ_{4} \sqrt{Δ p （ t ）} + θ_{5} 问_{}^{ˆ} （ t - 1 ）$ のパラメタ $θ_{3.}$ 、 $θ_{4}$ ，および $θ_{5}$ を推定します。この推定を実行する補助関数dynamicPipeEstを参照してください。

Th3 = 0 (3,1);Th4 = 0 (3,1);Th5 = 0 (3,1);[th3 (1) th4 (1) th5 (1)] = dynamicPipeEst (dp Q I1);% 3 4 5对区域1的估计[th3 (2), th4 (2), th5 (2)] = dynamicPipeEst (dp Q I2);% 3 4 5在区域2中的估计值[th3 (3), th4 (3), th5 (3)] = dynamicPipeEst (dp Q I3);% 3, 4, 5对区域3的估计

静的なポンプモデルの場合と異なり，動的な配管モデルは流量値への動的な依存を示しています。速度状態を変化させてモデルをシミュレートするため,控制系统工具箱の“LPV系统”ブロックを使用して仿真软件で区分的線形モデルが作成されます。金宝app金宝appSimulinkモデルLPV_pump_pipeと，このシミュレ，ションを実行する補助関数simulatePumpPipeModelを参照してください。

检查控制系统工具箱是否可用ControlsToolboxAvailable = ~isempty(ver(“控制”) &&许可证(“测试”，“Control_Toolbox”）;如果ControlsToolboxAvailable模拟动态管道模型。使用测量的压力值作为输入Ts = t(2)-t(1);Switch = ones(size(w));开关(I2) = 2;开关(I3) = 3;UseEstimatedP = 0;Qest_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);情节(t, Q t Qest_pipe)比较实测流量和预测流量其他的负载预保存的分段线性Simulink模型仿真结果金宝app负载DynamicOperationDataQest_pipeTs = t(2)-t(1);情节(t, Q t Qest_pipe)结束包含(“时间(s)”) ylabel (流量(Q) m^3/s)传说(“测量”，“估计”，“位置”，“最佳”)标题(“动态管道模型验证”）

C.動的なポンプ配管モデルの同定

動的なポンプ-配管モデルは，上記で特定したものと同じパラメタ( $θ_{3.} ， θ_{4} ， θ_{5}$ )を使用しますが,モデルのシミュレーションでは測定した値ではなく推定した圧力差を使用する必要があります。したがって，新しい同定の必要はありません。 $θ_{3.} ， θ_{4} ， θ_{5}$ の推定値によって，ポンプと配管のダ。

如果ControlsToolboxAvailable UseEstimatedP = 1;Qest_pump_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);情节(t, Q t Qest_pump_pipe)比较实测流量和预测流量其他的负载DynamicOperationDataQest_pump_pipe情节(t, Q t Qest_pump_pipe)结束包含(“时间(s)”) ylabel (流量Q (m^3/s))传说(“测量”，“估计”，“位置”，“最佳”)标题(“动态泵管模型验证”）

適合は，測定した圧力値を使って取得したものと事実上同一です。

D.動的な逆ポンプモデルの同定

測定したトルク値を前のトルクと速度の測定値に回帰させることにより，パラメ，タ， $θ_{6} ，．．．， θ_{10}$ を同様の方法で特定できます。しかし,結果の区分的線形モデルの完全な複数速度シミュレーションでは,データへの適切な適合が得られません。そのため，異なるブラックボックスのモデル化の方法を試します。これにはデタに当てはめる有理リグレッサをも非線形arxモデルの同定が含まれます。

在550个样本中使用前300个样本进行识别N = 350;sys3 = identifyNonlinearARXModel(Mmot,w,Q,Ts,N)

sys3 = 1输出2输入的非线性ARX模型输入:u1, u2输出:y1标准回归量对应阶数:na = [2] nb = [2 1] nk =[0 1]自定义回归量:u1(t-2)^2 u1(t)*u2(t-2) u2(t)^2非线性回归量:无模型输出在回归量中是线性的。采样时间:0.1秒状态:使用时域数据的NLARX估计。拟合估计数据:50.55%(模拟焦点)FPE: 1759, MSE: 3214

Mmot_est = sim(sys3，[w Q]);Mmot情节(t, t, Mmot_est)%比较测量和预测电机转矩包含(“时间(s)”) ylabel (“电机转矩(Nm)”)传说(“测量”，“估计”，“位置”，“最佳”)标题(“逆泵模型验证”）

残差の生成

モデルの"残差"を，測定信号と，モデルで生成された対応する出力との差として定義します。4

R1 = dp - dest;r2 = Q - Qest_pipe;r3 = Q - Qest_pump_pipe;

逆ポンプモデル残差の計算には,元の残差に大きな分散が見られるため,移動平均フィルターを使用してモデル出力に平滑化操作を適用します。

r4 = Mmot - movmean(Mmot_est，[1 5]);

学習残差のビュ:

图subplot(221) plot(t,r1) ylabel(“静态泵- r1”) subplot(222) plot(t,r2)动态管道- r2) subplot(223) plot(t,r3)动态泵管- r3)包含(“时间(s)”) subplot(224) subplot(t,r4)动态逆泵- r4)包含(“时间(s)”）

残差特徴の抽出

残差は，故障特定のための適切な特徴を抽出する信号です。利用できるパラメーター情報がないため,信号の最大振幅や分散などの信号プロパティのみから求められる特徴を考慮します。

PRBS入力実現を使用したポンプ-配管システムに対する20の実験セットを考えます。次の各モドにいて実験セットを繰り返します。

健全なポンプ
故障1:クリアランスの隙間にある摩損
故障2:@ @ンペラ，@ @出口にある少量の堆積物
故障3:@ @
故障4:@ @
故障5:ブレ，ド折損
故障6:キャビテ，ション
故障7:速度センサバアス
故障8:表示流量計バ
故障9:圧力センサバアス

実験デ，タの読み込み

url =“//www.tatmou.com/金宝appsupportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/MultiSpeedOperationData.mat”；websave (“MultiSpeedOperationData.mat”url);负载MultiSpeedOperationData生成操作模式标签标签= {“健康”，“ClearanceGapWear”，“ImpellerOutletDeposit”，.．.“ImpellerInletDeposit”，“AbrasiveWear”，“BrokenBlade”，“空化”，“SpeedSensorBias”，.．.“FlowmeterBias”，“PressureSensorBias”};

各アンサンブルおよび各操作モドにき残差を計算します。これには数分かかります。したがって、残差デ、タはデ、タファ、ルに格納されます。以下に示すように，helperComputeEnsembleResiduesを実行して残差を生成します。

% HealthyR = helpcomputeensembler余量(HealthyEnsemble,Ts,sys3,th1,th2,th3,th4,th5);%健康数据残差

从“helpercomputeensemblerremainder”运行加载预先保存的数据url =“//www.tatmou.com/金宝appsupportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/Residuals.mat”；websave (“Residuals.mat”url);负载残差

モ，ド区別の能力が最も高い残差の特徴は事前に知られていません。そのため,各残差について平均,最大振幅,分散,尖度,1ノルムなど,候補の特徴をいくつか生成します。すべての特徴が，健全なアンサンブルにおける値の範囲を使用してスケ，リングされます。

CandidateFeatures = {@mean, @ (x)马克斯(abs (x)), @kurtosis, @var, @ (x)和(abs (x))};featurename = {“的意思是”，“马克斯”，“峰度”，“方差”，“OneNorm”};%由收集的各故障模式残差生成特征表[HealthyFeature, MinMax] = helperGenerateFeatureTable(HealthyR, CandidateFeatures, featurename);Fault1Feature = helperGenerateFeatureTable(Fault1R, CandidateFeatures, featurename, MinMax);Fault2Feature = helperGenerateFeatureTable(Fault2R, CandidateFeatures, featurename, MinMax);Fault3Feature = helperGenerateFeatureTable(Fault3R, CandidateFeatures, featurename, MinMax);Fault4Feature = helperGenerateFeatureTable(Fault4R, CandidateFeatures, featurename, MinMax);Fault5Feature = helperGenerateFeatureTable(Fault5R, CandidateFeatures, featurename, MinMax);Fault6Feature = helperGenerateFeatureTable(Fault6R, CandidateFeatures, featurename, MinMax);Fault7Feature = helperGenerateFeatureTable(Fault7R, CandidateFeatures, featurename, MinMax);Fault8Feature = helperGenerateFeatureTable(Fault8R, CandidateFeatures, featurename, MinMax);Fault9Feature = helperGenerateFeatureTable(Fault9R, CandidateFeatures, featurename, MinMax);

各特徴テブルには20の特徴(残差信号ごとに5の特徴)があります。各テブルは，各実験から1ず，50(行)。

N = 50;每种模式下的实验数%FeatureTable = [.．.[HealthyFeature (1: N,:), repmat(标签(1)[N, 1]));.．.[Fault1Feature (1: N,:), repmat(标签(2),[N, 1]));.．.[Fault2Feature (1: N,:), repmat(标签(3),[N, 1]));.．.[Fault3Feature (1: N,:), repmat(标签(4),[N, 1]));.．.[Fault4Feature (1: N,:), repmat(标签(5),[N, 1]));.．.[Fault5Feature (1: N,:), repmat(标签(6),[N, 1]));.．.[Fault6Feature (1: N,:), repmat(标签(7),[N, 1]));.．.[Fault7Feature (1: N,:), repmat(标签(8),[N, 1]));.．.[Fault8Feature (1: N,:), repmat(标签(9),[N, 1]));.．.[Fault9Feature (1: N,:), repmat(标签(10),(N, 1]]];FeatureTable.Properties。VariableNames{结束}=“条件”；预览一些训练数据样本。disp(FeatureTable([2 13 37 49 61 62 73 85 102 120]，:))

Mean1非常刻薄的Mean3 Mean4 Max1 Max2 Max3 Max4 Kurtosis1 Kurtosis2 Kurtosis3 Kurtosis4 Variance1 Variance2 Variance3 Variance4 OneNorm1 OneNorm2 OneNorm3 OneNorm4条件  ________ _______ _______ _______ _______ ________ ________ _________ _________ _________ _________ __________ _________ _________ _________ _________ ________ ________ ________ ________ _________________________ 0.32049 0.6358 0.6358 0.74287 0.39187 0.53219 0.45647 0.6263 0.6263 0.53219 0.041795 0.18957 0.089492 0.089489 00.15642 0.56001 0.63541 0.63541 0.24258{'健康'}0.21975 0.19422 0.19422 0.83704 0 0.12761 0.12761 0.27644 0.19643 0.20856 0.033063 0.16761 0.16773 0.025119 0.19344 0.063167{'健康'}0.1847 0.25136 0.13929 0.1393 0.084162 0.72346 0.091488 0.091485 0.052552 0.063757 0.18371 0.18371 0.023845 0.10671 0.25065 0.25065 0.04848{'健康'}11 10 0.78284 0.94535 0.94535 0.75934 0.92332 11 0.80689 0.9978311 0.9574 {'Healthy'} -2.6545 0.90555 0.90555 0.86324 1.3037 0.7492 0.7492 0.97823 -0.055035 0.57019 0.57018 1.4412 1.4411 0.73128 0.90542 0.80573 3.2084 0.90585 1.0744 0.82197 0.30254 0.30254 -0.022325 0.21411 0.098504 0.0985 -0.010587 0.55959 0.29554 0.29554 0.066693 1.0432 0.3326 0.43326 0.13785 {'ClearanceGapWear'} -1.2149 0.53579 0.53579 1.0899 0.4739 0.47339 0.47339 0.29539 0.048137 0.17864 0.796480.48658 0.48658 0.25686 1.4066 0.5353 0.5353 0.26663 {'ClearanceGapWear'} -1.0949 0.44616 0.44616 1.143 0.56693 0.19719 0.19719 - 0.032603 0.32603 -0.0033592 0.43341 0.12857 0.12857 0.038392 0.2238 0.44574 0.093243 {'ClearanceGapWear'} -0.1844 0.16651 0.16651 0.87747 0.25597 0.080265 0.080265 0.49715 0.5019 0.29939 0.29938 0.5338 0.072397 0.058024 0.058025 0.1453 0.20263 0.16566 0.16566 0.14216 {'ImpellerOutletDeposit'} -3.0312 0.9786 0.75241 1.473 0.978390.97839 1.0343 0.0050647 0.4917 0.4917 0.89759 1.7529 0.9568 0.9568 0.8869 3.6536 0.97859 0.97859 0.62706 {'ImpellerOutletDeposit'}

分類器の設計

A.散布図を使用したモ，ドの分離可能性の可視化

特徴の視覚検査によって解析を始めます。そのためには，故障1:クリアランスの隙間にある摩損にいて考えます。この故障を検出するのに最も適した特徴を表示するため,“健康”と‘ClearanceGapWearのラベルを使って特徴の散布図を生成します。

T = FeatureTable(:，1:20);P = t变量;R = featutable . condition;I = strcmp(R，“健康”) | strcmp(R，“ClearanceGapWear”）;F =数字;gplotmatrix (P(我:),[],R (I)) f.Position (3:4) = f.Position (3:4) * 1.5;

明確には見えませんが，1列目と17列目の特徴で最もよく分離しています。これらの特徴をさらに詳しく解析します。

F =数字;名称= featutable . properties . variablenames;J = [1 17];流(间隙故障选择特征:%s\nstrjoin(名字(J),”、“）)

间隙故障选择特征:Mean1, OneNorm1

gplotmatrix (P (I, 17 [1]), [], R (I))

これでプロットには,特徴Mean1およびOneNorm1を使用して,クリアランスの隙間の故障モードから健全モードを分離できることが明確に示されます。同様の解析を各故障モドにいて実行できます。すべてのケースにおいて、故障モードを区別する特徴のセットを見つけることができます。したがって、故障動作の“検出”は常に可能です。ただし，同じ特徴に影響する故障タaaplプが複数あるため，故障の“特定”はより困難になります。たとえば,特徴Mean1 (r1の平均)と特徴OneNorm1 (r1の1ノルム)は多くの故障タイプで変化が見られます。同時に,センサーバイアスなど一部の故障は,多くの特徴において故障の特定が可能であり,比較的簡単に分離できます。

3。

图;I = strcmp(R，“健康”) | strcmp(R，“PressureSensorBias”) | strcmp(R，“SpeedSensorBias”) | strcmp(R，“FlowmeterBias”）;J = [1 4 6 16 20];选定特征指数百分比流(“传感器的选定特征”偏置:%s\n”strjoin(名字(J),”、“）)

传感器偏置的选择特征:Mean1, Mean4, Max2, Variance4, OneNorm4

gplotmatrix (P (I, J), [], R (I))

。故障の簡約セット(センサーバイアスのみ)に対し,特徴の簡約セットを使用して20メンバーのツリーバガー分類器に学習させます。

rng默认的再现率%Mdl = TreeBagger(20, FeatureTable(I，[21])，“条件”，.．.“OOBPrediction”，“上”，.．.“OOBPredictorImportance”，“上”）;图(oobError(Mdl)) xlabel(“树的数量”) ylabel (“误分类概率”）

誤分類のエラは3%未満です。したがって,故障の特定のサブセットを分類するため,さらに少ない数の特徴のセットを選択し,取り扱うことが可能です。

B.分類学習器アプリを使用した複数クラス分類器

前の節では散布図の手動検査に注目し，特定の故障タ▪▪▪プに▪▪▪いて特徴セットを小さくしました。この方法は手間がかかり，すべての故障タ。さらに自動化された方法ですべての故障モ，ドに対処できる分類器を設計することは可能でしょうか。统计和机器学习工具箱には数多くの分類器が用意されています。短時間でこれらの多くを試してその性能を比較する 1 つの方法として、分類学習器アプリを使用できます。

分類学習器アプリを起動して,ワークスペースから,新しいセッションの作業データとしてFeatureTableを選択します。ホルドアウト検証のためにデタの20%(各モドの10のサンプル)を確保しておきます。
メ@ @ンタブの[モデルタ电子邮箱プ]セクションで，[すべて]を選択します。次に[学習]ボタンを押します。
短時間で20の分類器の学習が行われます。これらの正確性は，履歴パネルでその名前の横に表示されます。線形SVM分類器の性能が最良であり,ホールドアウトサンプルでは86%の正確性が得られます。この分類器では“ClearanceGapWearの識別がやや難しく,40%の確率で“ImpellerOutletDepositとして分類されます。
性能のグラフィカルビューを得るには,メインタブの[プロット]セクションから混同行列プロットを開きます。プロットは，選択した分類器(ここでは線形SVM分類器)の性能を示します。

最も性能が良い分類器をワ，クスペ，スにエクスポ，トして，新しい測定値に対する予測に使用します。

まとめ

設計の優れた故障診断手法によって,サービスのダウンタイムとコンポーネントの交換費用を最小限に抑えることで運用コストを節約できます。この手法は,稼働するマシンのダイナミクスに関する適切な知識があると有利で,これとセンサー測定値を組み合わせて使用して,異なる種類の故障を検出し特定します。この例では，遠心ポンプの故障診断のための残差に基づく方法にいて説明しました。この方法は,モデル化タスクが複雑でモデルパラメーターに動作状態への依存が見られる場合に,パラメーターの推定と追跡に基づく方法の良い代替となります。

残差に基づく故障診断方法は以下の手順で行われます。

物理的な考慮事項またはブラックボックスのシステム同定手法を使用して,システムの測定可能な入力と出力の間のダイナミクスをモデル化します。
測定した信号とモデルで生成された信号の差として残差を計算します。必要に応じて，故障特定可能性を改善するために残差をさらにフィルタ，処理します。
各残差信号からピ，ク振幅，電力，尖度などの特徴を抽出します。
異常検出と分類の手法を用いて，故障の検出と分類のための特徴を使用します。
すべての残差および求められた特徴が，すべての故障の影響を受けるとは限りません。特徴ヒストグラムと散布図のビューによって,特定の故障タイプを検出するにはどの特徴が適しているかが明らかになります。特徴を選択し，故障特定におけるその性能を評価するこの過程は，反復手順となる場合もあります。

参考文献

Isermann,罗尔夫,故障诊断的应用。米odel-Based Condition Monitoring:Actuators, Drives, Machinery, Plants, Sensors, and Fault-tolerant System，第1版，斯普林格出版社柏林海德堡，2011年。

サポ，ト関数

静的なポンプ方程式のパラメ，タ，推定

函数[x1, x2, dest] = staticPumpEst(w, dp, I)静态泵参数估计在变化的速度设置% I:所选操作区域的样本指数。W1 = [0;w(我)];Dp1 = [0;dp (I)];R1 = [w1.]^ 2 w1);x = pinv(R1)*dp1;X1 = x(1);X2 = x(2);dest = R1(2:end，:)*x;结束

動的な配管のパラメ，タ，推定

函数[x3, x4, x5, Qest] = dynamicPipeEst(dp, Q, I)%dynamicPipeEst动态管道参数估计在变化的速度设置% I:所选操作区域的样本指数。Q = Q(i);dp = dp(I);R1 = [0;问(1:end-1)];R2 = dp;R2(R2<0) = 0;R2 =√(R2);R = [ones(size(R2))， R2, R1];去除不在状态的样品ii = find(I);J = find(diff(ii)~=1);R = R(2:end，:);R(j，:) = [];y = Q(2:结束);Y (j) = [];x = R\y;X3 = x(1);X4 = x(2);X5 = x(3); Qest = R*x;结束

LPV系统ブロックを使用した,ポンプ——配管モデルの動的な複数操作モードシミュレーション。

函数Qest = simulatePumpPipeModel(Ts,th3,th4,th5)动态管道系统分段线性建模。% Ts:采样时间% w:泵转速% th1, th2, th3为3种状态下的估计模型参数。此功能需要控制系统工具箱。ss1 = ss(th5(1)，th4(1)，th5(1)，th4(1)，Ts);ss2 = ss(th5(2)，th4(2)，th5(2)，th4(2)，Ts);ss3 = ss(th5(3)，th4(3)，th5(3)，th4(3)，Ts);抵消=排列([th3 (1) th3 (2), th3(3)]的[3 2 1]);OP = struct(“地区”，[1 2 3]');Sys = cat(3,ss1,ss2,ss3);sys。SamplingGrid = OP;assignin (“基地”，“sys”sys) assignin (“基地”，“抵消”，offset) MDL =“LPV_pump_pipe”；sim (mdl);Qest = logsout.get(“问”）;Qest = Qest. values;Qest = Qest. data;结束

逆ポンプダ@ @ナミクスの動的なモデルの特定。

函数syse = identifyNonlinearARXModel(Mmot,w,Q,Ts,N)识别2输入(w, Q)， 1输出(Mmot)数据的非线性ARX模型。%的输入:% w:转速% Q:流量% Mmot:电机转矩% N:要使用的数据样本数量%输出:syse:已识别的模型％这个函数使用系统识别工具箱中的NLARX估计器。Sys = idnlarx([2 2 1 0 1]，”，“CustomRegressors”,{“u1(2) ^ 2”，“u1 (t) * u2 (2)，“u2 (t) ^ 2”});data = iddata(Mmot，[w Q]，Ts);opt = nlarxOptions;opt.Focus =“模拟”；opt.SearchOptions.MaxIterations = 500;syse = nlarx(data(1:N)，sys,opt);结束