高速ダイナミクスの削減

高速ダイナミクス、すなわち高周波数のダイナミクスのソースを特定して減らすことで、モデルをリアルタイム対応にすることができます。リアルタイム対応モデルとは、リアルタイム マシン上でオーバーランを発生させずに許容可能な結果を出力するモデルです。この例では、線形化されたモデルの周波数応答と極配置を調べることによって高速ダイナミクスを特定する方法を説明します。また、高速ダイナミクスのソースを特定して除外する方法も示します。

高速ダイナミクスを削減する理由

高速ダイナミクスをもつモデルの計算コストは通常高くなります。高速ダイナミクスを除外すると、計算コストが下がり、固定ステップ、固定コストのシミュレーションに指定できる最小ステップ サイズが大きくなります。使用するステップ サイズが大きいと、モデルがリアルタイム対応となる可能性が高くなります。

周波数応答解析

"周波数応答" は、正弦波入力に対するシステムの定常状態応答を表します。線形システムの場合、正弦波入力は、周波数 ω が同じで振幅と位相 θ の異なる正弦波を出力します。

周波数解析は、与えられた周波数の範囲で振幅と位相がどのように変化するかを示します。小さな周波数の変化で振幅または位相が大きく変化するのは、システムに高速ダイナミクスがあることを示しています。この例では、振幅 (dB 単位の振幅) と位相が周波数の関数としてどのように変化するかを見ることができるボード線図を使用します。

極解析

高速極も、高速ダイナミクスがあることを示します。高速極は高速で応答または発振する極です。複素平面上で虚軸から左に離れた実数部をもつ極の応答速度は速くなります。実軸から離れた虚数部をもつ複素共役極対は高速で発振します。たとえば、-1500 の実極は -1000 の実極より応答速度が速く、-500 ± 1500i の複素共役極対は -500 ± 500i の複素共役極対より発振速度が速くなります。

状態空間モデルでは、極は A 行列の固有値です。この例では、状態空間モデルを決定してから A 行列値の固有値を計算しプロットすることによって極の速度を調べる方法を説明します。

モデルの線形化

この例のモデルは線形ではありません。周波数応答解析および極解析を実行する前に、モデルの平衡化 (つまり、線形化のための操作点を抽出して指定) と線形化を行います。

1. モデルを開いて調べます。MATLAB^® コマンド プロンプトに以下を入力します。

   %% Open the model
   open_system('ssc_hydraulic_actuator_digital_control')

   

   信号の生成、操作、経路指定、および可視化を行う各ブロックの他に、次のブロックがモデルに含まれています。

   • Controller — 制御システムの連続時間表現を定義する Transfer Fcn ブロック。モデルのコールバック関数は、伝達関数の分子 num と分母 den を、変数としてワークスペースに保存します。

   • — 連続時間制御システムの離散時間実装を展開する際の、計算遅延とサンプルアンドホールド関数に関連する遅延を表す Transport Delay ブロック。

   • Linearization I/O points — 平衡化の場合は閉ループ、線形化の場合は開ループとしてシステムを構成できるようにするサブシステム。ワークスペースで ClosedLoop を 1 に設定すると、モデルのコールバック関数はシステムを閉ループとして構成します。

   • Hydraulic Actuator — プラントの物理モデルを含むサブシステム。

2. システムの線形化に適した操作点を見つけます。モデルをシミュレートし、Simscape™ ログ ノードからデータを抽出したうえで、結果をプロットして調べます。

    モデル シミュレーションおよび結果のプロット用スクリプト

   %% Simulate the model and plot the simulation results
   % Simulate the model
   sim('ssc_hydraulic_actuator_digital_control')
   
   % Extract simulation results from the Simscape logging nodes
   simlog1 = simlog_ssc_hydraulic_actuator_digital_control;
   pCylA1 = simlog1.Hydraulic_Actuator.Hydraulic_Cylinder.Chamber_A.A.p.series;
   pCylB1 = simlog1.Hydraulic_Actuator.Hydraulic_Cylinder.Chamber_B.A.p.series;
   aValve1 = simlog1.Hydraulic_Actuator.Custom_2_Way_Valve.Calculate_Orifice_Area.PS_Saturation.O.series;
   
   % Plot simulation results 
   h3_ssc_hydraulic_actuator_digital_control=figure('Name',...
       'ssc_hydraulic_actuator_digital_control',...
       'Outerposition', [1244 673 576 513]);
   
   % Plot pressures
   ah(1) = subplot(2,1,1);
   plot(pCylA1.time,pCylA1.values,'LineWidth',6,'Color','b');
   hold on
   plot(pCylB1.time,pCylB1.values,'LineWidth',3,'Color','b');
   hold off
   grid on
   ylabel('Pressure (Pa)');
   title('Cylinder Pressures');
   legend({'Side A','Side B'},'Location','Best');
   
   % Plot valve activity
   ah(2) = subplot(2,1,2);
   plot(aValve1.time,aValve1.values,'LineWidth',3,'Color','b');
   grid on
   title('Valve Opening vs. Time');
   ylabel('Opening (m^2)');
   xlabel('Time (s)');
   linkaxes(ah,'x');
   
   % Create legend
   legend('Original','Location','southeast');
   

   

   カスタム双方向バルブは、シミュレーション時間 t が 2 ～ 3 秒のときに開いています。

3. モデルを平衡化します。バルブが開いている t = 2.5 秒の時間を操作点として使用し、閉ループ シミュレーションを実行します。

    モデル平衡化用スクリプト

   %% Trim the model
   assignin('base','ClosedLoop',1); % Close the feedback loop
   [t,x,y] = sim('ssc_hydraulic_actuator_digital_control');
   idx = find(t>2.5,1);                                            
   X = x(idx,:); 
   U = y(idx);
   

4. 関数 linmod を使用して連続時間モデルの開ループ構成を線形化し、状態変数 a、b、c、d をワークスペースに保存します。

    モデル線形化用スクリプト

   %% Linearize the model
   assignin('base','ClosedLoop',0); % Break the feedback loop
   [a,b,c,d]=linmod('ssc_hydraulic_actuator_digital_control',X,U);
   assignin('base','ClosedLoop',1); % Close the feedback loop
   

周波数応答解析と極速度解析の実行

1. ボード線図を生成します。

    ボード線図生成用スクリプト

   %% Calculate and plot the frequency response
   % Generate data
   npts = 100;
   w = logspace(-3,5,npts); 
   G = zeros(1,npts);
   for i=1:npts % Use negative feedback convention ( c:=-c, d:=-d)
       G(i) = (-c)*(2*pi*1i*w(i)*eye(size(a))-a)^-1*b +(-d);
   end
   
   % Create Bode plot figure
   h4_ssc_hydraulic_actuator_digital_control=figure...
       ('Name','hydraulic_actuator_digital_control');
   
   % Plot magnitude
   ah(1) = subplot(2,1,1);
   magline_h=semilogx(w,20*log10(abs(G)));
   xlim(ah(1),[w(1),w(end)]);
   grid on
   ylabel('Magnitude (dB)');
   title('Frequency Response (Hydraulic Actuator)');
   
   % Plot phase
   ah(2) = subplot(2,1,2);
   phsline_h=semilogx(w,180/pi*unwrap(angle(G)));
   set([magline_h,phsline_h],'LineWidth',1);
   ylabel('Phase (deg)'); 
   xlabel('Frequency (Hz)');
   grid on
   linkaxes(ah,'x');
   
   % Create legend
   legend('Original','Location','southwest');
   

   

   周波数 ω が 10² Hz から 10³ Hz になる間に、位相は約 600 度低下します。位相シフト θ の急激な変化は、システムに高速ダイナミクスがあることを示しています。

2. 関数 eig を使用して行列 a の固有値を計算し、複素平面に極をプロットします。

    行列 A の固有値計算およびプロット用スクリプト

   %% Calculate and plot the poles
   % Extract the poles
   poles = eig(a);
   
   % Extract the real and imaginary coordinates
   X1 = real(poles); 
   Y1 = imag(poles);
   
   % Create complex plane plot figure
   h5_ssc_hydraulic_actuator_digital_control = figure('Name',...
       'hydraulic_actuator_digital_control',...
       'OuterPosition',[668 150 1137 512]);
   
   % Create axes
   axes1 = axes('Parent',h5_ssc_hydraulic_actuator_digital_control,...
       'Position',[0.0315 0.112 0.932 0.815]);
   hold(axes1,'on');
   
   % Create title
   title('Complex Plane Plot');
   
   % Plot poles
   plot(X1,Y1,'DisplayName','Original Poles','Marker','*',...
       'MarkerSize',3,'LineStyle','none');
   
   % Limit and label axes
   xlim(axes1,[-3000 500]);
   xlabel({'Real','Axis'});
   ylim(axes1,[-2000 2000]);
   ylabel({'Imaginary',' Axis'});
   set(axes1,'XAxisLocation','origin','YAxisLocation','origin');
   box(axes1,'on');
   
   % Create legend
   legend(axes1,'show','Original','Location','southwest');
   

   

   発振している可能性のある 2 つの極の対を含め、6 つの高速極があります。

3. 極の対があることを確認します。関数 eigs を使用して、6 つの高速極の値をコマンド ウィンドウに出力します。

    極値出力用スクリプト

   %% Print fast-pole values
   eigs(a,6)
   

   ans =
   
      1.0e+03 *
   
     -2.0000 + 1.1547i
     -2.0000 - 1.1547i
     -0.4614 + 1.4208i
     -0.4614 - 1.4208i
     -1.0314 + 0.0000i
     -1.0000 + 0.0000i
   

   2 組の極の対があります。

高速ダイナミクスのソースの特定と排除

高速ダイナミクスの潜在的ソースについてモデルを調べます。

1. モデルを線形化するために、この例では関数 linmod を使用します。linmod のドキュメンテーションでは、Transport Delay ブロックを含むモデルの線形化にこの関数を使用しない方がよいとされています。Transport Delay ブロックのドキュメンテーションでは、線形化ルーチンのパデ近似によってモデルに動的状態が付加される可能性があることが指摘されています。このブロックがモデルの線形化をもたらす高速極のソースであるかどうかを判定します。

2. Transport Delay ブロックの影響を排除してモデルをシミュレートするには、このブロックをコメント スルーします。

    Transport Delay ブロックのコメント アウト用スクリプト

   %% Eliminate Transport Delay block from the network
   set_param('ssc_hydraulic_actuator_digital_control/Transport Delay',...
       'Commented','through')
   

   

   Transport Delay のアイコンの表示が薄くなり、このブロックがコメント スルーされたことがわかります。

3. Transport Delay ブロックに影響されていないモデルの周波数応答を調べるため、モデルの平衡化、線形化、およびシミュレーションを行い、その後ボード線図を更新します。

    モデルの平衡化と線形化、およびボード線図の更新用スクリプト

   %% Trim, Linearize, and update the Bode plot.
   % Trim
   assignin('base','ClosedLoop',1); % Close the feedback loop
   [t,x,y] = sim('ssc_hydraulic_actuator_digital_control');
   idx = find(t>2.5,1);                                            
   X = x(idx,:); 
   U = y(idx);
   
   % Linearize
   assignin('base','ClosedLoop',0); % Break the feedback loop
   [a,b,c,d]=linmod('ssc_hydraulic_actuator_digital_control',X,U);
   assignin('base','ClosedLoop',1); % Close the feedback loop
   
   % Generate data 
   npts = 100; w = logspace(-3,5,npts); G = zeros(1,npts);
   for i=1:npts % Use negative feedback convention ( c:=-c, d:=-d)
       G(i) = (-c)*(2*pi*1i*w(i)*eye(size(a))-a)^-1*b +(-d);
   end
   
   % Update the Bode plot figure
   figure(h4_ssc_hydraulic_actuator_digital_control);
   hold on
   ah(1) = subplot(2,1,1);
   hold on
   magline_h=semilogx(w,20*log10(abs(G)),'--');
   xlim(ah(1),[w(1),w(end)]);
   ah(2) = subplot(2,1,2);
   hold on
   phsline_h=semilogx(w,180/pi*unwrap(angle(G)),'--');
   set([magline_h,phsline_h],'LineWidth',1);
   linkaxes(ah,'x');
   legend('Original','No Transport Delay','Location','southwest');
   

   

   周波数 ω が 10² Hz から 10³ Hz になる間に、位相は約 250 度しか低下しません。

4. 高速極を計算してプロットします。

    行列 A の固有値計算およびプロット用スクリプト

   %% Calculate and plot the poles
   % Extract the poles
   poles = eig(a);
   
   % Extract the real and imaginary coordinates
   X1 = real(poles); 
   Y1 = imag(poles);
   
   % Update the complex plane plot figure
   figure(h5_ssc_hydraulic_actuator_digital_control)
   hold on
   plot(X1,Y1,'DisplayName','No Transport Delay','Marker','d',...
       'MarkerSize',5,'LineStyle','none');
   

   

   Transport Delay ブロックが、-2000 ± 1.1547i rad/s の発振極対がなくなっている原因です。

5. シミュレート結果をプロットして、元の結果と十分に一致しているかどうかを確認します。

    モデル シミュレーションおよび結果のプロット用スクリプト

   %% Plot the simulation results
   % Get simulation results
   simlog2 = simlog_ssc_hydraulic_actuator_digital_control;
   pCylA = simlog2.Hydraulic_Actuator.Hydraulic_Cylinder.Chamber_A.A.p.series;
   pCylB = simlog2.Hydraulic_Actuator.Hydraulic_Cylinder.Chamber_B.A.p.series;
   aValve = simlog2.Hydraulic_Actuator.Custom_2_Way_Valve.Calculate_Orifice_Area.PS_Saturation.O.series;
   
   % Update figure
   figure(h3_ssc_hydraulic_actuator_digital_control)
   hold on
   ah(1) = subplot(2,1,1);
   hold on
   plot(pCylA.time,pCylA.values,'LineWidth',3,'Color','r','LineStyle','-');
   hold on
   plot(pCylB.time,pCylB.values,'LineWidth',1,'Color','r','LineStyle','-');
   legend({'Side A','Side B'},'Location','Best');
   hold off
   ah(2)= subplot(2,1,2);
   hold on
   plot(aValve.time,aValve.values,'LineWidth',1,'Color','r','LineStyle','-');
   hold off
   legend({'Original','No Transport Delay'},'Location','southeast');
   

   

   結果は同じように見えます。

6. ズームして、より詳細に精度を評価します。

    ズーム イン用スクリプト

   %% Zoom the figure
   % Define the x-axis limits
   xStart = 0;
   xEnd = 10;
   xZoomStart1 = 1.78;
   xZoomEnd1 = 1.88;
   
   % Zoom in
   figure(h3_ssc_hydraulic_actuator_digital_control);
   ah(2);
   xlim([xZoomStart1, xZoomEnd1]);

   

   このレベルでは、修正モデルの結果にわずかな違いがみられます。しかしながら、シミュレーションは十分に正確であり、結果は実測データと理論データに基づく予測と一致しています。

7. モデルには、液圧圧縮性、すなわち油柱共鳴が含まれています。油柱共鳴が 2 番目の発振極対の原因となっていることを確認するため、圧縮性をオフにして、周波数応答解析と極解析を繰り返します。

    圧縮性排除、周波数応答および極解析用スクリプト

   %% Remove Compressibility
   model = 'ssc_hydraulic_actuator_digital_control';
   subsystem = 'Hydraulic Actuator';
   composite_block = 'Hydraulic Cylinder';
   
   CylA = [model,'/',subsystem,'/',composite_block,'/Chamber A'];
   CylB = [model,'/',subsystem,'/',composite_block,'/Chamber B'];
   set_param(CylA, 'Compressibility', '0')
   set_param(CylB, 'Compressibility', '0')
   
   %% Trim
   assignin('base','ClosedLoop',1); % Close the feedback loop
   [t,x,y] = sim('ssc_hydraulic_actuator_digital_control');
   idx = find(t>2.5,1);                                            
   X = x(idx,:); U = y(idx);
   
   %% Linearize
   assignin('base','ClosedLoop',0); % Break the feedback loop
   [a,b,c,d]=linmod('ssc_hydraulic_actuator_digital_control',X,U);
   assignin('base','ClosedLoop',1); % Close the feedback loop
   
   %% Update Bode plot
   % Generate data
   npts = 100; 
   w = logspace(-3,5,npts); 
   G = zeros(1,npts);
   for i=1:npts % Use negative feedback convention ( c:=-c, d:=-d)
       G(i) = (-c)*(2*pi*1i*w(i)*eye(size(a))-a)^-1*b +(-d);
   end
   
   figure(h4_ssc_hydraulic_actuator_digital_control);
   hold on
   ah(1) = subplot(2,1,1);
   hold on
   magline_h=semilogx(w,20*log10(abs(G)));
   xlim(ah(1),[w(1),w(end)]);
   % figure(h4_ssc_hydraulic_actuator_digital_control);
   % hold on
   
   ah(2) = subplot(2,1,2);
   hold on
   phsline_h=semilogx(w,180/pi*unwrap(angle(G)));
   set([magline_h,phsline_h],'LineWidth',1);
   linkaxes(ah,'x');
   legend('Original','No Transport Delay','No Compressibility',...
       'Location','southwest');
   
   %% Update complex plane plot
   poles = eig(a);
   
   figure(h5_ssc_hydraulic_actuator_digital_control)
   hold on
   X1 = real(poles); 
   Y1 = imag(poles);
   plot(X1,Y1,'DisplayName','No Compressibility','Marker','d',...
       'MarkerSize',8,'LineStyle','none');
   

   

   

   位相低下のさらなる減少は、高周波数のダイナミクスの縮小を反映しています。高速極は 2 つになっています。一方の高速極が虚軸からさらに離れて移動しても、高速ダイナミクスは減少しますが、これは、458.8 ± 1.4273i rad/sec で発振している極の対が除外されているためです。

8. 残りの 2 つの高速極の値をコマンド ウィンドウに出力します。

    極値出力用スクリプト

   %% Print the pole values
   eigs(a,2)
   

   two_fast_poles =
   
      1.0e+03 *
   
      -2.6767
      -1.0000
   

   -2677 rad/s の高速極は、油圧サブシステム内の 2 つのオリフィス コンポーネントによって生じる負荷質量と油圧の減衰に対応しています。これらのダイナミクスは、シミュレーション結果の中核をなす部分です。-1000 rad/s の高速極は、コントローラーの分母 0.001 s+1 に対応します。伝達関数はコントローラーの設計に不可欠です。重要なシステムレベルの動作を変えずに、これ以上動的モードを排除することはできません。

9. シミュレート結果をプロットして、元の結果と十分に一致しているかどうかを確認します。

    シミュレーション結果のプロット用スクリプト

   %% Print simulation results
   % Get simulation results
   simlog3 = simlog_ssc_hydraulic_actuator_digital_control;
   pCylA = simlog3.Hydraulic_Actuator.Hydraulic_Cylinder.Chamber_A.A.p.series;
   pCylB = simlog3.Hydraulic_Actuator.Hydraulic_Cylinder.Chamber_B.A.p.series;
   aValve = simlog3.Hydraulic_Actuator.Custom_2_Way_Valve.Calculate_Orifice_Area.PS_Saturation.O.series;
   
   % Update results
   figure(h3_ssc_hydraulic_actuator_digital_control)
   hold on
   ah(1)= subplot(2,1,1);
   hold on
   plot(pCylA.time,pCylA.values,'LineWidth',1,'Color','y','LineStyle','--');
   hold on
   plot(pCylB.time,pCylB.values,'LineWidth',1,'Color','y','LineStyle','--');
   legend({'Side A','Side B'},'Location','Best');
   hold off
   ah(2)= subplot(2,1,2);
   hold on
   plot(aValve.time,aValve.values,'LineWidth',1,'LineStyle','--');
   legend({'Original','No Transport Delay','Compressibility'},'Location','southeast');
   hold off
   
   % Zoom out
   figure(h3_ssc_hydraulic_actuator_digital_control);
   ah(2);
   xlim([xStart, xEnd]);
   

   

   更新されたモデルの精度は許容可能であるように見えます。

10. ズームして、より詳細に精度を評価します。

     ズーム イン用スクリプト

    %% Zoom the figure
    % Zoom in
    figure(h3_ssc_hydraulic_actuator_digital_control);
    ah(2);
    xlim([xZoomStart1, xZoomEnd1]);
    

    

    このレベルでは、修正モデルの結果にわずかな違いしか加わっていないことがわかります。シミュレーションの精度は許容可能です。

参考

ブロック

• Transport Delay

関数

• eig | eigs | linmod

関連する例

• ステップ サイズの決定
• 計算コストの推定
• 計算コストの削減
• 数値的剛性の軽減
• ゼロクロッシングの削減

詳細

• モデルを準備する目的
• リアルタイム モデルの準備のワークフロー
• モデルの線形化