【MATLAB】摩擦下の制御系の相軌道を描いてみた

はじめに

この記事では、非線形摩擦によって乱されてしまった制御系の挙動がどんなものになるのか、相軌道という切り口で図示してみたいと思います。
最終的には、MATLABでこんな絵が描けます。

準備

本題に入る前に、問題設定と表現の紹介をしておきます。

題材

図のように平面に置かれた物体の、位置制御問題を題材にします。
運動方向は１次元で右向きを正にとり、物体の座標を $x[m]$ 、初期位置は $x_0$、目標位置は $x=0$ 、物体の質量を $m[kg]$ とします。
物体は平面から摩擦力 $f[N]$ を受け、我々はこの物体に制御入力として力 $F[N]$ を与えられるものとします。

この系の運動方程式は以下で与えられます。

m\ddot{x} = F+f

摩擦モデル

物体に作用する摩擦は、以下のモデルで表される非線形摩擦のみとします。
これはクーロン摩擦（静止摩擦と動摩擦）のモデルとして一般に用いられているものです。
なお、ここでは簡単化のために最大静止摩擦と動摩擦の大きさを共に $f_0[N]$ としています。
速度 $\dot{x} = 0$ の時は物体を動かす力 $F$ を相殺しようとする力（上限 $f_0$ ）が、$\dot{x} \neq 0$ の時は物体が減速する方向に一定の力 $f_0$ が作用します。

f(t) = \begin{cases}-F(t) & (\dot{x}=0, |F|\leq f_0)\\
-sgn(F(t))f_0 & (\dot{x}=0, |F|> f_0)\\
-sgn(\dot{x}(t))f_0 & (\dot{x}\neq 0)
\end{cases}

制御入力

制御入力 $F$ は、以下に示すPID制御（$K_p, K_i,K_d$：定数ゲイン、$\eta$：積分量）で与えるものとします。

\begin{align}
F(t) &= -(K_p x(t) + K_i \eta(t) + K_d \dot{x}(t))\\
\eta(t) &= \int_0^t x(\tau) d\tau
\end{align}

相空間と相軌道

以上で定義されたシステムは微分方程式で記述され、その挙動は解である時間関数 $x(t)$、および $\dot{x}(t), \eta(t)$ の変化の過程として表現できます。
そこで、この３つの関数で張った三次元空間を考えると、解の時間変化は初期値から始まる軌跡を描きます。
この三次元空間が相空間、解の描く軌道が相軌道に相当します。

本題

上記のシステムのシミュレーションをMATLABで実装します。

シミュレーション

ode45を用いて、複数の初期値に対するシミュレーションを回してプロットします。
なお、現象の比較のために摩擦の有無を切り替える変数fswを導入しています。

main.m

%% 変数の初期化とfigureのclose
clear;
close all;

%% 定数の設定
global Ki Kp Kd m f0 fsw%グローバル変数として定義
fsw = 1;% 0 or 1

% wall 3
Ki = 0.5;
Kp = 1;
Kd = 1;
m = 1;
f0 = 0.05;

%% シミュレーション
x0 = -0.5:0.02:0.5;
figure;
for i=1:length(x0)
  [t, x] = ode45('odeFcn', [0:0.1:30], [0;x0(i);0]);
  hold on
  if x0(i)>0
      plot3(x(:,2), x(:,3), x(:,1),'r');
  else
      plot3(x(:,2), x(:,3), x(:,1),'b');
  end
end
grid on
xlabel('x')
ylabel('dx')
zlabel('\eta')

%% single plot
[t, x] = ode45('odeFcn', [0:0.1:50], [0;0.1;0]);
figure;
subplot(1,2,1)
plot(t,x(:,2));
grid on;
xlabel('t');
ylabel('x');

subplot(1,2,2)
plot3(x(:,2), x(:,3), x(:,1));
grid on
xlabel('x')
ylabel('dx')
zlabel('\eta')

odeFcn.m

function dx = odeFcn(t, x)
    global Ki Kp Kd m f0 fsw

    F = -( Kp*x(2) + Ki*x(1) + Kd*x(3) );
    f = - sign(x(3))*f0 - (1-sign(x(3))^2) * (sign(F) * min(abs(F),f0));

    dx  = [x(2);
           x(3);
           (F + fsw*f)/m];
end

現象１：仮に摩擦なしとした場合

非線形摩擦がない場合には、適当なゲインによって（値によっては振動を伴いながら）$x=0$ へと指数的に収束させられることが知られています。
上記のコードで変数 fsw=0 とした場合の応答と、ある初期値から得られる相軌道の例を図に示します。
見慣れた線形システムの応答ですね。

一連の初期値から得られる相軌道を重ねて表示すると、次のようになります。
安定性を感じますね。

現象２：摩擦が存在する場合

一方、クーロン摩擦の影響下では、見慣れない応答波形が得られる場合があります。
具体的には次の図のように、途中で運動を停止する期間を含む応答です。

何が起きているのでしょうか？

相軌道の方を見ていると、停止している($\dot{x}=0$)間にも状態は変化していて、z軸の積分量 $\eta$ のみが増加していることが分かります。
偏差が小さくP制御成分だけではクーロン摩擦を超えるだけの力が出せない状態になると、このような停滞が発生するということです。
加えてこの例においては、行き過ぎが発生してから運動の方向を反転する際に、クーロン摩擦の符号が逆転するため、より顕著に影響が現れます。

また、一連の相軌道を確認してみると、$\dot{x} = 0$ の平面上に断層が形成されていることが分かります。
このようにして、非線形摩擦による相軌道の幾何形状の変化を確認することが出来ました。

まとめ

クーロン摩擦の影響下では、シンプルな制御対象とPID制御においてさえも、複雑な応答になってしまう事例を紹介しました。
筆者自身、制御工学に触れて間もない頃に実機でこの現象に出くわして、困惑した記憶があります。
そんな現象も相軌道として描画することで、視覚的にある種の特徴を捉えることが出来ました。
考察の切り口として、何かの参考になれば幸いです。