[PYTHON] Introduction of data-driven controller design method

Introduction

Feedback control has been used to "control" everything from home appliances to industrial equipment. Those that perform this control are called controllers, and various design methods have been considered. This time, I would like to introduce a relatively new approach called "data-driven controller design method".

In addition, I will leave the theoretical part to the dissertation and write mainly on using it. So, I would like to explain it with the code of MATLAB and Python.

Up to data-driven control

Roughly speaking, there are three approaches to controller design:

1. Manual parameter adjustment (limit sensitivity method, etc.)
2. Model-based controller design (optimal regulator, etc.)
3. Others ($ H_ {\ infinty} $ control, data-driven control, etc.)

Since 1 is manual, it is not difficult to imagine that the design is difficult because it requires trial and error. Therefore, as in 2, a method was proposed in which the * mathematical model to be controlled was exactly obtained * and the model was used for designing. The task of finding this model exactly is called "identification." By describing the optimization problem using the identified model, we are now able to obtain the theoretically optimal controller. However, it is the musician who wants this ** model exactly **, and it takes a lot of time and effort to reconcile theory and experiment. Therefore, other methods were considered based on the premise that ** models cannot be obtained exactly **. Data-driven control is one of them.

What is data-driven control?

According to the website of Professor Kaneko of the University of Electro-Communications, what is data-driven control?

A method of designing a controller directly from the data and behavior of the device to be controlled when it is driven.

is. In other words, the controller can be designed without creating a model, as long as there is signal data during driving.

There are currently two well-known data-driven controller design methods [^ excuse IFT].

[^ excuseIFT]: Actually, historically, non-proof control and IFT need to be explained, but I will omit them.

This time, I wrote the code about VRFT because the optimization problem can be solved easily.

Sample code execution environment

• MATLAB 2017b
  • Control System Toolbox
  • System Identification Toolbox
• Python 3.7.5 (Anaconda environment)
  • matplotlib
  • numpy
  • scipy
  • python-control

Note that the Python code in the explanation omits import statements. See [Whole Code](#Whole Code) for working code.

VRFT procedure

There are only three things a designer needs to do to use VRFT.

1. Determination of design specifications
2. Acquisition of input / output data
3. Solve the optimization problem according to the algorithm

This time, we will explain the controller design of the speed control system using a motor as an example.

Determining design specifications

First, let's decide the design specifications. The following three design specifications can be determined.

For example, if you want to realize the response of the time constant $ \ tau = 0.02 $ [s] by the PID controller,

\begin{align}
M(s) &= \frac{1}{0.02 s + 1} \\&= \frac{50}{s + 50}\\
\boldsymbol{\beta}(s) &= \begin{bmatrix}1& \frac{1}{s}& s\end{bmatrix}^\mathrm{T}
\end{align}

Set as follows. Also, considering that there is a high possibility that noise will be mixed in the data in the high frequency range, it is often sufficient to set a low-pass filter for the frequency weight.

The above code is as follows.

vrft.m

%Sampling time and operators
Ts = 0.001;
s = tf('s');        %Laplace operator s
z = tf('z', Ts);    %Time advance operator z

%Reference model M
tau = 0.02;
M = c2d(1/(tau*s + 1), Ts); %First-order lag system discretized by zero-order hold

%Weight function
gW = 100;
W = c2d(gW/(s + gW), Ts);   %First-order lag system discretized by zero-order hold

%Control structure beta
beta = minreal([1; Ts/(1 - z^-1); (1 - z^-1)/Ts]); %PID controller

vrft.py

# ===Determining design specifications===
#Sampling time
Ts = 0.001

#Reference model M
tau = 0.02
M = ctl.c2d(ctl.tf([1], [tau, 1]), Ts)  #First-order lag system discretized by zero-order hold

#Weight function W
gW = 100
W = ctl.c2d(ctl.tf([gW], [1, gW]), Ts)  #First-order lag system discretized by zero-order hold

#Control structure beta
A, B, C, D = abcd_normalize(
    [[1., 0], [0, 0]],
    [[1.], [-1.]],
    [[0, 0], [Ts, 0], [0, -1/Ts]],
    [[1.], [0], [1/Ts]])
beta = ctl.ss(A, B, C, D, Ts)
del A, B, C, D

Acquisition of input / output data

Next, the input / output data $ \ {u_0 (t), y_0 (t) } $ is acquired from the control target. In the case of VRFT, the data is acquired so that the following conditions are met.

• Be sure to get the open loop response (only $ P (s) $ without control)
• If possible, use the white input signal $ u_0 (t) $

As the whiteness signal, it is recommended to use the M-sequence signal that is often used for system identification. In the case of MATLAB, M-sequence signals can be easily generated using the idinput function of the System Identification Toolbox. In the case of Python, it can be generated by Scipy's max_len_seq function. However, it seems that signals of a different series from MATLAB are generated, so in this Python version, I made my own function prbs () and used it.

Here, a model of the controlled object (motor in this case) is required to perform the simulation [^ note1]. So, I originally thought that it was necessary to explain the model of the motor, but Manao's Great article I found items / ed383a30bbe0c7b9534a). Therefore, we will omit modeling here and assume that the motor can be represented by a primary system.

[^ note1]: Note that this is only needed for data acquisition, no model is used in the design. (Actually, it can be designed as clearvars P after data acquisition.)

Therefore, in the simulation, data can be acquired with the following code.

vrft.m

%Input signal used for data acquisition u(m-sequence signal)
n = 15;                                     %Number of stages
T = 2^n - 1;                                %Number of data per cycle
p = 15;                                     %Number of data cycles
N = T*p;                                    %The number of data
u0 = idinput([T 1 p],'prbs',[0,1],[-1,1]);   %Signal vector

%Input signal power spectral density phi_u
phi_u = 1;                  %Input signal is assumed to be white

%Controlled model P
Tp = 0.74;
Kp = 1.02;
P = c2d(Kp/(Tp*s + 1), Ts);

%Output signal y0
y0 = lsim(P, u0);

vrft.py

# ===Acquisition of input / output data===
def prbs(n, p):
    # matlab compatible PRBS signal generator

    taps = {3: [0, -1], 4: [0, -1], 5: [1, -1], 6: [0, -1], 7: [0, -1],
            8: [0, 1, 6, -1], 9: [3, -1], 10: [2, -1], 11: [8, -1],
            12: [5, 7, 10, -1], 13: [8, 9, 11, -1], 14: [3, 7, 12, -1],
            15: [13, -1], 16: [3, 12, 14, -1], 17: [13, -1], 18: [10, -1]}
    N = (2**n - 1)*p
    x = np.ones(n, dtype=np.int8)
    u = np.zeros(N, dtype=np.int8)
    tap = taps[n]
    for i in range(N):
        u[i] = x[-1]
        x0 = x[tap[0]] ^ x[tap[1]]
        x = np.roll(x, 1)
        x[0] = x0
    return u


#Input signal used for data acquisition u(m-sequence signal)
n = 15
T = 2**n - 1
p = 15
N = T*p
# u0, _ = max_len_seq(n, length=N)  #The signal series is different from MATLAB
u0 = prbs(n, p)
u0 = -(2.*u0 - 1.)

#Input signal power spectral density phi_u
phi_u = 1.

#Controlled model P
Tp = 0.74
Kp = 1.02
P = ctl.c2d(ctl.tf([Kp], [Tp, 1]), Ts)

#Output signal y0
y0, t0, _ = ctl.lsim(P, u0)

Solve the optimization calculation according to the algorithm

Finally, the signal is generated according to the VRFT algorithm and the optimization calculation is solved.

This is a minimization problem for the merit function $ J_ \ mathrm {MR} (\ boldsymbol {\ rho}) $, because VRFT aims to respond as close to the reference model $ M (z) $ as possible.

\begin{align}
J_\mathrm{MR}(\boldsymbol{\rho}) &= \left\| \left( \frac{P(z) C(z, \boldsymbol{\rho})}{1 + P(z) C(z, \boldsymbol{\rho})} - M(z) \right) W(z) \right\|_2 ^2
\end{align}

Here, $ \ boldsymbol {\ rho} $ represents the control parameter, $ C (z, \ boldsymbol {\ rho}) = \ boldsymbol {\ rho} ^ {\ rm T} \ boldsymbol {\ beta } (z) There is a $ relationship. However, since this formula uses the controlled model $ P (z) $, it needs to be rewritten with I / O data. The theoretical explanation is cut off, but the problem of minimizing the evaluation function $ J_ \ mathrm {MR} $ is the following, and the problem of minimizing $ J_ \ mathrm {VR} ^ N $ is asymptotic. Is equivalent to.

\begin{align}
J_{\rm VR}^{N}(\boldsymbol{\rho}) &= \frac{1}{N} \sum_{t=1}^{N}(u_L(k) - \boldsymbol{\rho}^{\rm T} \boldsymbol{\varphi}(t))^2
\end{align}

Here, the time series data $ u_L (t) $ and $ \ boldsymbol {\ rho} $ can be generated by the following formula using the filter $ L (z) $.

\begin{align}
u_L(t) &= L(z) u_0(t)\\
y_L(t) &= L(z) y_0(t)\\
%\boldsymbol{\varphi}(k) &= \boldsymbol{\beta}(z) (\tilde{r}_L(t) - y_L(t))
e_L(t) &= (M^{-1} - 1)y_L(t)\\
\boldsymbol{\varphi}(k) &= \boldsymbol{\beta}(z) e_L(t)
\end{align}

However, the filter $ L (z) $ shall satisfy the following.

\begin{align}
|L(z)|^2 &= \frac{|1-M(z)|^2 |M(z)|^2 |W(z)|^2}{\phi_u(\omega)} &\forall \omega
\end{align}

Then, $ J_ \ mathrm {VR} ^ N (\ boldsymbol {\ rho}) $ does not use the controlled model $ P (z) $, and the input / output data $ \ {u_0 (t), y_0 (t) It can be described only with) } $. Also, since this equation is linear with respect to the control parameter $ \ boldsymbol {\ rho} $, it can be solved by the least squares method.

From the above, the part that solves the optimization problem is the following code.

vrft.m

%Prefilter L
L = minreal(M*(1 - M)/phi_u);

%Input signal passed through filter ul
ul = lsim(L, u0);

%Pseudo-error signal el
el = lsim(L*(M^(-1) - 1), y0);

%Control output before parameters phi
phi = lsim(beta, el);

%Optimal parameters rho
rho = phi\ul;               %Solve the least squares method in matrix form(mldivide)

%Designed controller C
C = minreal(rho.' * beta);  %Ask for a controller

%Evaluation function Jmr
Jmr = mean(ul - phi * rho); %Check the evaluation function in matrix format

vrft.py

# ===Designed by VRFT===
#Prefilter L
L = ctl.minreal(M*(1 - M)/phi_u)

#Input signal passed through filter ul
ul, _, _ = ctl.lsim(L, u0)

#Pseudo-error signal el
el, _, _ = ctl.lsim(ctl.minreal(L*(M**-1 - 1)), y0.flatten())

#Control output before parameters
phi, _, _ = ctl.lsim(beta, el.flatten())

#Optimal parameters rho
solution = np.linalg.lstsq(phi, ul, rcond=None)
rho = solution[0]

#Designed controller C
C = ctl.ss([0], [0, 0, 0], [0], rho.T, Ts) * beta   #Ask for a controller

#Evaluation function Jmr
Jmr = np.mean(ul - np.dot(phi, rho))    #Check the evaluation function in matrix format

Confirmation of performance

With the above procedure, it seems that the controller could be designed using only the input / output data without using the model to be controlled. The optimal control parameters $ \ rho $ obtained are as follows. You can see that MATLAB and Python give almost the same results.

[^ ideal]: $ M (z) = \ frac {C ^ * (z) P (z)} {1 + C ^ * (z) P (z)} A control that satisfies $ $ C ^ * $ Value (this is called an ideal controller)

Now, let's check if the response of the reference model is really realized by this control. The code to check is as follows.

vrft.m

%Whole system with controller G
G = minreal(feedback(P*C, 1));

%Step response
fig1 = figure('name', 'Step plot');
stepplot(G, M);
legend({'$$\frac{CP}{1+CP}$$', '$M$'},...
    'Interpreter', 'latex', 'FontSize', 14, 'Location', 'southeast');

%Bode plot display
fig2 = figure('name', 'Bode plot of controller');
bodeplot(G, M, {1,100});

vrft.py

# ===Confirmation of performance===
#Whole system with controller G
G = ctl.minreal(ctl.feedback(P*C, 1))

#Step response
plt.figure()
plt.title("Step response of closed loop system")
timeRange = np.arange(0, 0.12 + Ts, Ts)
ym, _ = ctl.step(M, timeRange)
yg, _ = ctl.step(G, timeRange)
plt.plot(timeRange, ym, timeRange, yg)
plt.xlabel("Time [s]", usetex=True)
plt.ylabel("Velocity [V]")
plt.legend(['Reference model', 'Closed loop system'], loc='lower right')

#Bode plot display
plt.figure()
plt.title("Bode plot of closed loop system")
ctl.bode(M, G)
plt.legend(['Reference model', 'Closed loop system'], loc='lower left')
plt.show()

Here, only the step response is shown as an image. You can see that the response of the reference model can be reproduced with indistinguishable accuracy.

Whole code

Finally, I will show you the whole code I have introduced.

%VRFT design script
%Designed for the speed control system of the motor.

% Copyright (c) 2019 larking95(https://qiita.com/larking95)
% Released under the MIT Licence 
% https://opensource.org/licenses/mit-license.php

%%Initialization
clearvars;
close all;

%%Determining design specifications
%Sampling time and operators
Ts = 0.001;
s = tf('s');        %Laplace operator s
z = tf('z', Ts);    %Time advance operator z

%Reference model M
tau = 0.02;
M = c2d(1/(tau*s + 1), Ts); %First-order lag system discretized by zero-order hold

%Weight function
gW = 100;
W = c2d(gW/(s + gW), Ts);   %First-order lag system discretized by zero-order hold

%Control structure beta
beta = minreal([1; Ts/(1 - z^-1); (1 - z^-1)/Ts]); %PID controller

%%Acquisition of input / output data
%Input signal used for data acquisition u(m-sequence signal)
n = 15;                                     %Number of stages
T = 2^n - 1;                                %Number of data per cycle
p = 15;                                     %Number of data cycles
N = T*p;                                    %The number of data
u0 = idinput([T 1 p],'prbs',[0,1],[-1,1]);   %Signal vector

%Input signal power spectral density phi_u
phi_u = 1;                  %Input signal is assumed to be white

%Controlled model P
Tp = 0.74;
Kp = 1.02;
P = c2d(Kp/(Tp*s + 1), Ts);

%Output signal y0
y0 = lsim(P, u0);

%%Designed by VRFT
%Prefilter L
L = minreal(M*(1 - M)/phi_u);

%Input signal passed through filter ul
ul = lsim(L, u0);

%Pseudo-error signal el= 
el = lsim(L*(M^(-1) - 1), y0);

%Control output before parameters phi
phi = lsim(beta, el);

%Optimal parameters rho
rho = phi\ul;               %Solve the least squares method in matrix form(mldivide)

%Designed controller C
C = minreal(rho.' * beta);  %Ask for a controller

%Evaluation function Jmr
Jmr = mean(ul - phi * rho); %Check the evaluation function in matrix format

%%Confirmation of performance
%Whole system with controller G
G = minreal(feedback(P*C, 1));

%Step response
fig1 = figure('name', 'Step plot');
stepplot(G, M);

%Bode plot display
fig2 = figure('name', 'Bode plot of controller');
bodeplot(G, M, {1,100});

# -*- coding: utf-8 -*-
"""
VRFT design script
Designed for the speed control system of the motor.

Copyright (c) 2019 larking95(https://qiita.com/larking95)
Released under the MIT Licence 
https://opensource.org/licenses/mit-license.php
"""

import matplotlib.pyplot as plt
import numpy as np
import control.matlab as ctl

# from scipy.signal import max_len_seq
from scipy.signal import abcd_normalize


def prbs(n, p):
    # matlab compatible PRBS signal generator

    taps = {3: [0, -1], 4: [0, -1], 5: [1, -1], 6: [0, -1], 7: [0, -1],
            8: [0, 1, 6, -1], 9: [3, -1], 10: [2, -1], 11: [8, -1],
            12: [5, 7, 10, -1], 13: [8, 9, 11, -1], 14: [3, 7, 12, -1],
            15: [13, -1], 16: [3, 12, 14, -1], 17: [13, -1], 18: [10, -1]}
    N = (2**n - 1)*p
    x = np.ones(n, dtype=np.int8)
    u = np.zeros(N, dtype=np.int8)
    tap = taps[n]
    for i in range(N):
        u[i] = x[-1]
        x0 = x[tap[0]] ^ x[tap[1]]
        x = np.roll(x, 1)
        x[0] = x0
    return u


# ===Determining design specifications===
#Sampling time
Ts = 0.001

#Reference model M
tau = 0.02
M = ctl.c2d(ctl.tf([1], [tau, 1]), Ts)  #First-order lag system discretized by zero-order hold

#Weight function W
gW = 100
W = ctl.c2d(ctl.tf([gW], [1, gW]), Ts)  #First-order lag system discretized by zero-order hold

#Control structure beta
A, B, C, D = abcd_normalize(
    [[1., 0], [0, 0]],
    [[1.], [-1.]],
    [[0, 0], [Ts, 0], [0, -1/Ts]],
    [[1.], [0], [1/Ts]])
beta = ctl.ss(A, B, C, D, Ts)
del A, B, C, D

# ===Acquisition of input / output data===
#Input signal used for data acquisition u(m-sequence signal)
n = 15
T = 2**n - 1
p = 15
N = T*p
# u0, _ = max_len_seq(n, length=N)  #The signal series is different from MATLAB
u0 = prbs(n, p)
u0 = -(2.*u0 - 1.)

#Input signal power spectral density phi_u
phi_u = 1.

#Controlled model P
Tp = 0.74
Kp = 1.02
P = ctl.c2d(ctl.tf([Kp], [Tp, 1]), Ts)

#Output signal y0
y0, t0, _ = ctl.lsim(P, u0)

# ===Designed by VRFT===
#Prefilter L
L = ctl.minreal(M*(1 - M)/phi_u)

#Input signal passed through filter ul
ul, _, _ = ctl.lsim(L, u0)

#Pseudo-error signal el
el, _, _ = ctl.lsim(ctl.minreal(L*(M**-1 - 1)), y0.flatten())

#Control output before parameters
phi, _, _ = ctl.lsim(beta, el.flatten())

#Optimal parameters rho
solution = np.linalg.lstsq(phi, ul, rcond=None)
rho = solution[0]

#Designed controller C
C = ctl.ss([0], [0, 0, 0], [0], rho.T, Ts) * beta   #Ask for a controller

#Evaluation function Jmr
Jmr = np.mean(ul - np.dot(phi, rho))    #Check the evaluation function in matrix format

# ===Confirmation of performance===
#Whole system with controller G
G = ctl.minreal(ctl.feedback(P*C, 1))

#Step response
plt.figure()
plt.title("Step response of closed loop system")
timeRange = np.arange(0, 0.12 + Ts, Ts)
ym, _ = ctl.step(M, timeRange)
yg, _ = ctl.step(G, timeRange)
plt.plot(timeRange, ym, timeRange, yg)
plt.xlabel("Time [s]", usetex=True)
plt.ylabel("Velocity [V]")
plt.legend(['Reference model', 'Closed loop system'], loc='lower right')

#Bode plot display
plt.figure()
plt.title("Bode plot of closed loop system")
ctl.bode(M, G)
plt.legend(['Reference model', 'Closed loop system'], loc='lower left')
plt.show()

in conclusion

This time, while omitting the theory completely, I introduced the outline of data-driven controller design and the actual code for one of them, VRFT. Control Engineering I don't know as much as the other participants in the Advent Calendar, but I hope it helps everyone to be interested in data-driven control.

Also, if you have any concerns, please feel free to comment.

References

1. Introduction to VRFT (Campi and Savaresi's site)
2. Introduction of FRIT (Introduction article by Professor Kaneko)
3. [Basics of system identification](https://www.amazon.co.jp/%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0%E5%90% 8C% E5% AE% 9A% E3% 81% AE% E5% 9F% BA% E7% A4% 8E-% E8% B6% B3% E7% AB% 8B-% E4% BF% AE% E4% B8% 80 / dp / 4501114800) (Book by Dr. Adachi)
4. System identification of DC motor using Arduino and MATLAB / Simulink (Manao )
5. NumPy for Matlab users(Scipy.org)