[PYTHON] [Introduction to PID] I tried to control and play ♬

PID control seems to be a technology that has been around for a long time. So, I carefully formulated it and played with it, so I will summarize it.

【reference】 ①PID controller ② [PID control](https://ja.wikipedia.org/wiki/PID%E5%88%B6%E5%BE%A1#:~:text=PID%E5%88%B6%E5%BE%A1 % EF% BC% 88% E3% 83% 94% E3% 83% BC% E3% 82% A2% E3% 82% A4% E3% 83% 87% E3% 82% A3% E3% 83% BC% E3 % 81% 9B% E3% 81% 84,% E6% 96% B9% E6% B3% 95% E3% 81% AE% E3% 81% 93% E3% 81% A8% E3% 81% A7% E3% 81% 82% E3% 82% 8B% E3% 80% 82)

It can be easily represented by the following graph. From the above reference ① By Arturo Urquizo - Here, each quantity is defined as follows.

\begin{align}
r(t)&;Target amount\\
y(t)&;After operation, the current amount of target plants and processes\\
e(t)&=r(t)-y(t);Current amount-Residual calculated with target amount\\
u(t)&=P+I+D;Operation amount\\
          &;Operations on target plants and processes\\
　　      &;The value of the current that flows to raise and lower the temperature and the angle for operating the handle and valve. ..\\
P&=K_pe(t);Proportional control\\
I&=K_i\int_0^t e(\tau )d\tau;Integral control\\
D&=K_d\frac{de(t)}{dt};Derivative control\\
\end{align}

Here, in the case of computer control, $ y (t) $ is a discrete value, so it is better to define integral control and differential control by numerical integration and numerical differentiation as follows. At sampling time $ t = t_n $, place as follows. The numerical differentiation may be of a higher order, but the following may be used in consideration of practicality and reactivity.

\int_0^te(\tau)d\tau = \Sigma_{i=0}^ne(t_i)\Delta t\\
(\frac{de(t)}{dt})_{t_n}=\frac{e(t_n)-e(t_{n-1})}{\Delta t}\\
\Delta t = t_n-t_{n-1} ; const=Set to 1

Try to drop it in the code

First, the control targets are as follows. In other words, it is an image of temperature control of electric furnaces.

def ufunc(M,u=0,t=0):
    M0=0.2
    m=0.05
    return M - m*(M-M0-u)

The controller side considers the following PID control. Part of the code is based on the following. 【reference】 Try PID control with Python

t = 100
params_list = [(0,0,0)]  #Kp,Ki,Set Kd
for kp,ki,kd in params_list:
    Kp,Ki,Kd=kp,ki,kd
    M0 = 2
    M=M0
    goal = 1.00
    e = 0.00 
    e1 = 0.00 
    u=0
    es=0
    x_list = []
    y_list = []
    u_list = []
    x_list.append(0)
    y_list.append(M)
    u_list.append(u)
    for i in range(1,t,1):
        M1 = M
        e1 = e
        e = goal - y_list[i-1] #Deviation (e)=Objective value (goal)-Present value
        es += e #Integral is addition
        u = Kp * e + Ki * es + Kd * (e-e1) #Derivative is subtraction
        M = ufunc(M,u,i) #Response from the target
        x_list.append(i)
        y_list.append(M)
        u_list.append(u)
    plt.hlines([goal], 0, t, "red", linestyles='dashed') #The goal is displayed with a red dotted line
    plt.plot(x_list, y_list, color="b") #The time change of the response is displayed in blue
    plt.title('Kp={}, ki={}, kd={}'.format(Kp, Ki, Kd))
    #plt.plot(x_list, u_list, color="black") #Operation amount is displayed in black
    plt.ylim(-0.5, goal*2+0.2) #Graph height adjustment
    plt.pause(0.2)
    plt.savefig("output/Kp,Ki,Kd/imageKp{}Ki{}Kd{}M0{}.png ".format(Kp,Ki,Kd,M0))    
    plt.close()

See the target reaction with params_list = [(0,0,0)]

When t = 0, if M = 2 and room temperature M0 = 0.2, it gradually approaches M = 0.2 as follows.

This system is expressed by the following equation.

\frac{dM}{dt}=-m(M-M0)

The solution is

M=M0+(2-M0)e^{-mt}\\
t=When 0, M=2

Ziegra Nichols limit sensitivity method

See the critical vibration of the object According to Wikipedia, if Ki = Kp / Ti and Kd = Kp · Td, the optimum parameters are determined by the marginal sensitivity method of Ziegra Nichols as follows.

Here, Ku and Tu are defined as follows. "The proportional gain Kp is gradually increased from 0, and when the control amount reaches the stable limit and the constant amplitude vibration is maintained, the increase of Kp is stopped. The proportional gain at this time is Ku (limit sensitivity). , If the vibration cycle is Tu, When the critical vibration is calculated for the above object, it is calculated as Tu = 2.0 when Ku = 39 as shown below.

If you rewrite the above table with Ki, Kd,

result

In fact, if you try with the above parameters, P control Kp19.5Ki0.0Kd0.0M02dt1 PI control Kp17.55Ki10.57Kd0.0M02dt1 PID control Kp23.4Ki23.4Kd5.85M02dt1 In other words, PID control fails for some reason This seems to be due to the rounding error.

dt precision extension and recalculation

So, extend dt so that it can be calculated in arbitrary increments as shown below.

def ufunc(M,u=0,t=0,dt=1):
    M0=0.2
    m=0.05
    return M - m*(M-M0-u)*dt

t = 100
dt=0.1

Ku=39
Tu=2
params_list =[(0.5*Ku,0,0),(0.45*Ku,0.54*Ku/Tu,0),(0.6*Ku,1.2*Ku/Tu,3*Ku*Tu/40)]
for kp,ki,kd in params_list:
    Kp,Ki,Kd=kp,ki,kd
    ...
    for i in range(1,int(t/dt),1):
        ...
        e = goal - y_list[i-1] #Deviation (e)=Objective value (goal)-Last realized value
        es += e*dt
        u = Kp * e + Ki * es + Kd * ((e-e1)/dt )
        M = ufunc(M,u,i,dt)

        x_list.append(i*dt)
        y_list.append(M)
        u_list.append(u)
        ...
    plt.savefig("output/Kp,Ki,Kd,dt/imageKp{}Ki{}Kd{}M0{}dt{}.png ".format(Kp,Ki,Kd,M0,dt))    
    plt.close()

Re-executed with dt = 0.1, I was able to draw well as expected. PID control Kp23.4Ki23.4Kd5.85M02dt0.1 Since the big change ends at t = 0-10, it is now possible to calculate the dt step in detail, so let's recalculate with the step dt = 0.01. P control Kp19.5Ki0.0Kd0.0M02dt0.01 PI control Kp17.55Ki10.57Kd0.0M02dt0.01 PID control Kp23.4Ki23.4Kd5.85M02dt0.01

Graph reproduction of PID control

Try to reproduce the following graph of English Wikipedia. It seems that the conditions are not drawn, so when I draw variously under the most likely conditions, the Ki dependency becomes almost similar under the following conditions.

• If you narrow down the conditions more, you can reproduce it completely, but I stopped here. In this graph, you can see the properties when the Ki parameter is changed sufficiently. The conditions are as follows The control target is The value was set to m = 0.4, which is larger than the above, and M0 = 0, which corresponds to room temperature.

def ufunc(M,u=0,t=0,dt=1):
    M0=0
    m=0.4
    return M + m*u*dt - m*(M-M0)*dt

The calculation is done with the same parameter values as shown. Three graphs were drawn at the same time with a legend. The drawing range of the graph is adjusted.

t = 20
dt=0.01
params_list =[(1,0.5,1),(1,1,1),(1,2,1)]
for kp,ki,kd in params_list:
    Kp,Ki,Kd=kp,ki,kd
    M0 = 0
    M=M0
    goal = 1.00
    e = 0.00 
    e1 = 0.00 
    u=0
    es=0
    x_list = []
    y_list = []
    u_list = []
    x_list.append(0)
    y_list.append(M)
    u_list.append(u)
    for i in range(1,int(t/dt),1):
        M1 = M
        e1 = e
        e = goal - y_list[i-1] #Deviation (e)=Objective value (goal)-Last realized value
        es += e*dt
        u = Kp * e + Ki * es + Kd * ((e-e1)/dt )
        M = ufunc(M,u,i,dt)
        x_list.append(i*dt)
        y_list.append(M)
        u_list.append(u)
    plt.hlines([goal], 0, t, "red", linestyles='dashed') #The goal is displayed with a red dotted line
    plt.plot(x_list, y_list,label="Kp{}Ki{}Kd{}".format(Kp,Ki,Kd)) 
    plt.legend()
    plt.ylim(-0.1, 1.5) #Adjust the height of the graph
    plt.pause(0.2)
plt.savefig("output/Kp,Ki,Kd,dt/imageKp{}Ki{}Kd{}M0{}dt{}2.png ".format(Kp,Ki,Kd,M0,dt))    
plt.close()

On the other hand, the Kp dependency can be drawn as follows, and it cannot be reproduced without selecting the parameters. And the Kd dependency seems to be almost reproduced as follows.

About the control target

This time, the control target was set by imagining the temperature of the electric furnace. In actual control, you need to control something more dynamic. In that case, I think that it can be realized by rewriting ufunc () introduced this time. Also, depending on the control target, I think that there are targets that cannot be PID controlled. We need to pursue further the limits of that area.

Summary

・ I tried playing with PID control ・ Control parameters were determined by the marginal sensitivity method of Ziegra-Nichols. ・ I tried to reproduce the graph posted on Wikipedia ・ I was able to explicitly define a function that describes the controlled object. ・ Introduced sampling time dt

・ I would like to actually apply it to temperature control and motor control.