# pid
# 系統
$$
d \frac{d}{d t} x{\left (t \right )} + k x{\left (t \right )} + m \frac{d^{2}}{d t^{2}} x{\left (t \right )} = F_{c}
$$
其解
$$
x{\left (t \right )} = C_{1} e^{\frac{t \left(- d - \sqrt{d^{2} - 4 k m}\right)}{2 m}} + C_{2} e^{\frac{t \left(- d + \sqrt{d^{2} - 4 k m}\right)}{2 m}} + \frac{F_{c}}{k}
$$

由上式，所以我們可以知道，如果沒有任何外力，system 會回到原點(-exponential)。
其實這樣看係數很難看出什麼，最好在laplace domain 做
或是直接試試，調動每個參數的影響

import sympy as sp

t, kp, ki, kd, p, i, d, fc= sp.symbols('t, k_p, k_i, k_d, p, i, d, F_c')
m,k,d= sp.symbols('m,k,d')
x=sp.Function('x')
e=sp.Function('e')(t) # e=x_target-x

#%% get sol
eq= sp.Eq(m*x(t).diff(t,2)+d*x(t).diff(t)+k*x(t), fc)
sol= ( sp.dsolve(eq,x(t)) )
print(sp.simplify(sol))

#%% simulate behavior
#x(0)=-10, x'(0)=0
sol= sp.dsolve( eq.subs(fc,0), x(t), ics={x(0):-10, x(t).diff(t).subs(t,0):0} ) #for no input

sol= sol.subs([(m,1),(k,10),(d,1)]) #change parameter as you want
sp.plot(sol.rhs, xlim=(0,5), ylim=(-10,10)) #x

大致上是 k大反應快，震動大。 d大反應慢，震動小

# with control
接下來 考慮controller
force ($F_c$) 是我們的control input。我們要算出到底要施多少力才能達到我們希望的效果
if using pid controller，則系統變成

$$
e=x_{target}-x
\\d \frac{d}{d t} x{\left (t \right )} + k x{\left (t \right )} + m \frac{d^{2}}{d t^{2}} x{\left (t \right )} = k_{d} \frac{d}{d t} e{\left (t \right )} + k_{i} \int e{\left (t \right )}\, dt + k_{p} e{\left (t \right )}
$$

這個是加了controller的系統

$x_{target}$ 當然可以是時變的(function of t)，但這裡假設我們很簡單的希望$x$回到原點，所以

for $x_{target}=0,\ e=0-x$
$$
d \frac{d}{d t} x{\left (t \right )} + k x{\left (t \right )} + m \frac{d^{2}}{d t^{2}} x{\left (t \right )} = k_{d} \frac{d}{d t} \left(- x{\left (t \right )}\right) + k_{i} \int \left(- x{\left (t \right )}\right)\, dt - k_{p} x{\left (t \right )}
$$
整理後
$$
m\frac{d^2 x(t)}{dt^2}+(d+k_d)\frac{dx(t)}{dt}+(k+k_p)x(t)+k_{i} \int x{\left (t \right )}\, dt=0
$$
可以看到，$k_p, k_d$ 其實就是分別把 等效彈力係數($k'=k+k_p$)和 等效阻尼係數 加大。
所以我們可以透過施加適當的外力來改變系統動態，可以讓e.g.讓系統反應加快（by $k_p$）