一条悬浮架(m)和上方车体(1/2M)刚性连接，因为分析的是半边悬浮架，所以车体负载取一半。
悬浮架上有4个线圈，通过输入电流产生向上的悬浮力$F_{1-4}$。线圈的位置用L14确定,用向量L_v代表L14组成的向量，L的值代表线圈和中心的偏移量，考虑到左右则L_v=[1.05 0.35 -0.35 -1.05];。悬浮架上两个传感器，位置和1、4号线圈重合。传感器的位移$z_1,z_2$取向上为正方向，取系统的状态变量为$[z_1, z_2, \stackrel{\cdot}{z_1},\stackrel{\cdot}{z_2}]$，初始状态为[-0.02 -0.02 0 0], 平衡位置为[-0.01, -0.01, 0, 0]。
因为$J\frac{\mathrm{d} \omega}{\mathrm{d} t}=M$，$\omega$取逆时针为正。$\omega=(\stackrel{\cdot}{z_2}-\stackrel{\cdot}{z_1})/2L$，则
$$
\begin{equation}
J\frac{\mathrm{d} \omega}{\mathrm{d} x} = \frac{J}{2L}(\stackrel{\cdot\cdot}{z_2}-\stackrel{\cdot\cdot}{z_1})
\tag{2}
\end{equation}
$$
力矩为:
$$
\begin{equation}
M = -\sum^{4}F_iL_i
\tag{3}
\end{equation}
$$
所以
$$
\begin{equation}
\frac{J}{2L}(\stackrel{\cdot\cdot}{z_2}-\stackrel{\cdot\cdot}{z_1})=-\sum^{4}F_iL_i
\tag{4}
\end{equation}
$$
用MATLAB求解$(1)、(4)$的联立，这里记$\frac{J}{2L}$为一个新常数$A_0$，求解结果为：
% sF 和 sFL分别代表两个求和式
syms x1 x2 m M sF sFL g A0
eq1 = 0.5*(0.5*M+m)*(x1+x2) == sF - (0.5*M+m)*g;
eq2 = A0*(x2-x1) == -sFL;
[ddz1,ddz2] = solve([eq1 eq2],[x1 x2])
% 输出：
ddz1 =
(4*A0*sF + M*sFL + 2*m*sFL - 2*A0*M*g - 4*A0*g*m)/(2*A0*(M + 2*m))
ddz2 =
-(M*sFL - 4*A0*sF + 2*m*sFL + 2*A0*M*g + 4*A0*g*m)/(2*A0*(M + 2*m))
所以得到系统的状态方程为：
$$
\left{\begin{matrix}
\stackrel{\cdot }{x_1} = x_3\
\stackrel{\cdot }{x_2} = x_4\
\stackrel{\cdot }{x_3} = \frac{1}{A_0M}(2A_0\sum F + m\sum FL - A_0Mg - 2A_0mg) \
\stackrel{\cdot }{x_4} = -\frac{1}{A_0M}(-2A_0\sum F + m\sum FL + A_0Mg + 2A_0mg)
\end{matrix}\right.
\tag{5}
$$
其中力是电流、线圈距离和角度的函数，可以写成$F_i=f_i(z_1,z_2,i)$，具体公式从前人处直接继承Simulink的相关子系统。
根据状态空间搭建Simulink框图：
将该系统在平衡位置线性化，使用Simulink 线性化工具，得到的结果为： $$ A = \begin{pmatrix} 0& 0& 1& 0\ 0& 0& 0& 1\ 5.5\times 10^5& -5.5\times 10^5& 0& 0\ -5.5\times 10^5& 5.5\times 10^5& 0&0 \end{pmatrix} \ B = \begin{pmatrix} 0& 0\ 0 & 0\ 187.9& -187.4\ -187.4& 187.9 \end{pmatrix} $$ 该结果保存在linsys_equib.mat中。
基于该线性模型涉及LQR控制器
>>K = lqr(linsys1,eye(4),1)
K =
1.0e+03 *
6.4569 0.5870 0.0902 0.0845
0.5870 6.4569 0.0845 0.0902
此时没有添加参考值，所有状态变量最终趋近于0。
对$\vec{x}$添加参考值$\vec{r}$，即令$\vec{u}=-K(\vec{x}-\vec{r})$。在这里$\vec{r}=[1,1,0, 0]^T\bullet r$，原系统和LQR控制器形成了新的以$r$值为输入，$\vec{x}$为输出的1输入4输出模型。其状态空间为： $$ {\vec{\stackrel{\cdot}{x}}} =A_{cl}\vec{x}+BK\begin{pmatrix} 1\1\0\0 \end{pmatrix}r $$ 分析它的阶跃响应，首先将其转化为传递函数形式。
input:
Acl = linsys1.A-linsys1.B*K;
BB = linsys1.B*K*[1 1 0 0]';
[b,a] = ss2tf(Acl,BB,eye(4),zeros(4,1))
output:
b =
1.0e+09 *
0 0 0.0000 0.0069 3.5826
0 0 0.0000 0.0069 3.5826
0 0.0000 0.0069 3.5826 0.0000
0 0.0000 0.0069 3.5826 -0.0000
a =
1.0e+09 *
0.0000 0.0000 0.0013 0.0923 1.7913
所以对$x_1$来说，$X_1(s)/R(s)=\frac{2次多项式}{4次多项式}$的形式。求阶跃响应时，将$R(s)=\frac{1}{s}$乘过去，用MATLAB的redisue函数求它的留数（r）、极点（p）和余项（k）。代表着可以把$X_1(s)$拆分为如下形式： $$ X_1(s)=\frac{a}{s}+\sum \frac{r_i}{s-p_i}+k(s) $$ 在对上式做拉普拉斯逆变换后，a就是系统的稳态值。 程序如下： input:
num = b;
den = [a,0];
[r,p,k]=residue(num, den)
step(step(b(1,:),a))
output:
r =
-0.0000
0.0000
173.7282
-175.7282
2.0000
p =
1.0e+03 *
-1.2538
-0.8785
-0.0406
-0.0401
0
k =
[]
上式结果说明在极点是0的时候对应的留数是2，也就是说系统的稳态响应是输入量的2倍，想获得-0.01的稳态值，就应该给-0.005的参考值作为输入。
X_1的阶跃响应
