已知线性方程：dq^2/dt=3306.63q-33.7i+0.1p, di/dt=-623.34i+53.19u, q是位置，i是电流，p是干扰，y是在有噪声v的情况下的位置测量值，y=q+v，求开环传递函数，并画出系统方框图，写出matlab求的过程

首先将dq^2/dt转化为dq/dt乘以q的形式，有：

dq/dt = (3306.63q - 33.7i - 0.1p) / q

将di/dt代入上式中得到：

dq/dt = (3306.63q - 33.7i - 0.1p) / q
di/dt = -623.34i + 53.19u

将y=qv带入得到：

y' = dq/dt * v = (3306.63q - 33.7i - 0.1p) / q * v

将i代入上式中，得到：

y' = (3306.63q - 33.7(-y'/v) - 0.1p) / q * v

整理可得开环传递函数：

G(s) = y(s) / p(s) = (3306.63s) / (s^2 + 33.7s + 0.1(3306.63))

系统方框图如下：

+-----------------+
p(s) | |
-------->| G(s) |-----> y(s)
| |
+-----------------+

Matlab代码如下：

syms s q i p v y;
eq1 = diff(q,t) == (3306.63*q - 33.7*i - 0.1*p) / q;
eq2 = diff(i,t) == -623.34*i + 53.19*u;
eq3 = y == q*v;
eq4 = subs(eq1,i,-y/v);
eq5 = subs(eq4,eq3);
G = simplify(y/p);
G_open = subs(G,p,1);
G_open = subs(G_open,v,1);
G_open = subs(G_open,u,0);
G_open = subs(G_open,y,1);
G_open = collect(G_open,q);
G_open = simplify(G_open);
num = coeffs(G_open,q,'All');
den = [1 33.7 0.1*3306.63];
G_open = tf(num,den);
step(G_open);