2018 22nd International Conference on System Theory, Control and Computing (ICSTCC)

State-feedback control algorithms for a CNC

machine
Dora Sabau, Petru Dobra
Automation Department
Technical University of Cluj Napoca
Cluj-Napoca, Romania

Abstract—The controller of a small size Computer Numerical In this paper, a MIMO model with two inputs and two
Controlled (CNC) system is designed and tested. The system is outputs is taken into consideration. In chapter 2, based on
considered to have a Multi Input Multi Output(MIMO) model the acquired data consisting on the input and the output
and two control design techniques are presented for it: Pole-
placement and Linear Quadratic Regulator(LQR). The LQR signals from the machine, a mathematical model is obtained.
problem is proposed to be solved by choosing the weights matrices The parameter identification results, obtained using System
using the total energy of the system. The simulation results are Identification ToolBox from Matlab, are presented in chapter
evaluated and the controller that meets the desired behavior is 3. Chapter 4 focuses on designing and simulating two state-
implemented on the real machine. feedback control techniques: pole-placement method and Lin-
Index Terms—LQR, Pole-placement, energy, MIMO, CNC,
control ear Quadratic Regulator(LQR) and chapter 5 and 6 show the
experimental results and the conclusions.
I. I NTRODUCTION
Numerical control machines started to be used since 1950s.
Because of the rapid development of the industry, these ma- II. MATHEMATICAL MODEL FOR CNC
chines were improved annually. Currently, Computer Numeri-
cal Control(CNC) machines are very popular in the industrial
field in domains such as: aeronautical industry, energy and The workstation in the laboratory, which represents a small-
power technology and manufacturing industry. These machines size CNC machine is shown in Fig.1.
provide a high level of automation and consistent motion
control, accuracy, high precision and power, fast machining in
accordance with modern tools [1]. CNC systems are, regularly,
compound of: the machine itself and the command system.
The control system incorporates new techniques which reduce
the allocated time to a project and allow a complex control
of the process. The requirements these control systems have
nowadays are: precise control of the position, with safe veloc-
ity control and good acceleration and deceleration properties
[2]. In order to meet these requirements, the controller that is
used has to give good precision and to have high resolution
[3], [4].
As shown in literature [5], [6], conventional Proportional-
Integrative-Derivative(PID) controller, State-Feedback con-
trol, Feed-Forward control, Linear Quadratic Regulator(LQR), Fig. 1: CNC machine
Fuzzy Logic Control or even more advanced control tech-
niques such as optimized PD control using genetic algorithms
or adaptive control, are used to control CNC machines. The positioning system compounds of three axes of trans-
In order to obtain the control law for positioning systems lation: X, Y and Z. Z-axis is driven by a stepper motor
with translation axes, it is important to have a mathematical and has the role to place the tool at the right distance on
model for the process [7]. For this, two particular directions the vertical direction. X and Y-axis are driven by BrushLess
are taken into consideration. The first one is to consider the Direct Current(BLDC) motors which receive commands from
system as a decoupled one, which means each axis is handled MC206X controller. The controller is connected to the PC via
as an independent process. The second approach takes into an USB link, as shown in Fig.2 and can be programmed using
consideration the coupled-axis, so the process is treated in TrioBasic programming Language, which gives the possibility
terms of Multi-Input-Multi-Output(MIMO) model [8]. to develop control commands fast and easy.

978-1-5386-4444-7/18/$31.00 ©2018 IEEE 615

 into consideration only the time constant Tm, representing the
mechanical part of the drive:
KMx KMy
Hωx ux = ; Hωy uy = (1)
TMx s + 1 TMy s + 1
The second transfer function shows the relationship between
the speed and the position and is represented by an integrator:
1 1
Hθx ωx = ; Hθy ωy = ; (2)
s s
The mutual effect between the two axes can be approxi-
mated with proportional gains:
Fig. 2: Controller connections Hωy ux = Kωy ux ; Hωx uy = Kωx uy (3)

The purpose of this paper is to consider a mathematical The behavior analysis and the control design algorithm will
model for 2D positioning system(taking into consideration X be made using state space model:
and Y-axis) in order to design a control law for it. (
ẋ = Ax + Bu
(4)
y = Cx + Du
The notations used are as follows:
• x-state vector
• A-the state matrix
• B-the input matrix
• u-the vector of the input signals
• C-the output matrix
• y-the vector of the output signals
• D-the direct transfer
The state space variables are considered the measured
variables: the rotor speed and position for each axis. The input
signals vector is represented by the command signal ’u’ for
Fig. 3: Block diagram representation for the model used for each axis. The outputs are considered only the variables that
parameters estimation need to be controlled: the position of each axis.
   
The system can be represented, as shown in Fig.3, by a x1 (t) ωx (t)
Multi-Input-Multi-Output process with two inputs and two x2 (t)  θx (t) 
x(t) = x3 (t) = ωy (t)
   (5)
outputs, as it follows:
• Inputs: ux , uy -represent the drive command, where index x4 (t) θy (t)
x or y suggest the corresponding axis.
 
ux (t)
• Outputs: θx , θy -represent the position of rotors of the u(t) = (6)
uy (t)
motors that drive each axis.  
Also, we have the designations that suggest the type of the θ (t)
y(t) = x (7)
signal or the relationship between signals: θy (t)
• ωx , ωy : angular velocity; Transforming the transfer function model represented in
• Hωx ux , Hωy uy , Hθx ωx ,Hθy ωy the transfer functions from Fig.3 into a state space model with the previous notations and
u to ω and from ω to θ, for each axis considerations, we obtain:
• Hωx uy , Hωy ux the mutual effect between the two axes
1
− TM 0 0 0
 
The block diagram represented in Fig.3 illustrates the x

simplified scheme that represents the dependency between  1 0 0 0

A= 1
 (8)
electrical and mechanical part of the system. The whole system  0 0 − TM 0
y
can be represented by two coupled subsystems, each one 0 0 1 0
corresponding to one axis. The mutual effect between the axes  KM 
x
Kωx uy
appears because of the mechanical connection between the  M
T x
two.  0 0 
B= (9)

KMy 
The transfer function from the command u to the speed ω, Kωy ux T My

for each axis, can be approximated by a first order term, taking 0 0

616
    
0
1 0 0 −40.85 0 0 0
C= (10)
0
0 0 1  1 0 0 0
A= ;
   0 0 −87.81 0
0 0 0 0 1 0
D= (11)
0 0  
1.05 0.267
III. PARAMETER I DENTIFICATION  0 0 
B = 103 ∗ 
0.0245 2.188 ;
 (17)
Using TrioBasic commands, we have the measured signal
obtained by following a path at the maximum speed of the 0 0
 
rotor of each axis. The acquired data signals are represented 0 1 0 0
C= ;
in Fig.4, with the mention that the signals are scaled for a 0 0 0 1
good representation of all the signals in the same figure.
 
0 0
D= ;
0 0
IV. C ONTROLLER DESIGN
A. Pole Placement
Closed-loop pole assignment is a common technique used to
design controllers for MIMO processes. Specifying the closed
loop eigenvalues, a state-feedback matrix K is computed, so
that the closed loop system satisfies the condition [11]:

(A − BK)x = λi Ix (18)

where λi are the desired eigenvalues.

Also, to eliminate the stationary heading error when tracking
a reference signal, a pre-filter matrix F is calculated using the
following equation:
Fig. 4: Acquired data for each X and Y-axis
F = −(C(A − BK)−1 B)−1 (19)
Using the models described above and applying the iden-
tification methods [9] based on minimization of prediction To allocate the eigenvalues [-80 -80 -120 -120], the matrices
error and the least squares algorithm, the model parameters K and F have the following values:
are obtained using System Identification Toolbox in Matlab  
[10]. 0.1510 9.1107 −0.0013 −0.1110
K= (20)
Once the fitting, autocorrelation and cross-correlation tests −0.0017 −0.1019 0.513 4.388
are passed, the estimated model is obtained and transformed
into transfer functions using ’Zero Order Hold’ discretisation
 
9.1107 −0.1110
method. So, we have: F = (21)
−0.1019 4.388
• For X-axis:
1054 The step response is simulated and represented in Fig.5.
Hωx ux = (12)
s + 40.85
1
Hθx ωx = (13)
s
• For Y-axis:
2188
Hωy uy = (14)
s + 87.81
1
Hθy ωy = (15)
s
• Mutual effect:
Hωy ux = Kωy ux = 24.46; Hωx uy = Kωx uy = 26.65
(16)
The state space model will have: Fig. 5: Closed-loop step response using pole-placement

617
B. LQR Problem +2xT P T (Ax + B(−R−1 B T P x)) = 0 (31)
Another way to design a controller for a MIMO process is which gives the Ricatti equation mentioned above.
using LQR problem. This problem is based on optimal control, Since we have the Ricatti equation, we have solved the
which is concerned with operating a system at minimum cost. optimal control problem. The feedback control law is:
In this case, the system dynamics is described using linear
1
differential equations (state space model) and the cost function u = − R−1 B T (2P (t)x) = −R−1 B T P (t)x = −Kx (32)
is represented by a quadratic function as the following one: 2
Z ∞ where K is the gain matrix.
Like in the pole-placement case, if we want to implement a
J= (xT (t)Qx(t) + uT (t)Ru(t))dt (22)
0 tracking controller, it is necessary to compute a matrix N the
where Q defines the weights on the state and is a semi-defined reference signals is multiplied by, so that the stationary error
symmetric matrix and R defines the weights on the control is null.
input.
The aim is to minimize the cost function J and to find the N = −(C(A − BK)−1 B)−1 (33)
optimal control law u* that should be used [12]. In many cases, the weight matrix Q or R is chosen through
Based on the scientific articles in the literature, it turns out successive attempts, so that the results fulfill the requirements.
that regardless of the matrix Q and matrix R, the minimum This paper propose a method to chose the weights of matrix
cost function is obtained by solving a Ricatti equation [13]: Q using the energy of the system [8], [14].
Having:
−Ṗ (t) = Q − P T BR−1 B T P + AT P + P A (23) • the general differential equation of a motor, correspond-

This equation is obtained by solving the Hamilton-Jacobi- ing to each axis:

Bellman (HJB) equation (23) with the boundary condition in ∂2θ ∂θ
(24). J 2
+D + Kθ = u (34)
∂t ∂t
∂J ∗ (x(t), t) • the energy of the subsystem, consisting of each axis
+ minu H = 0 (24) separately:
∂t

J ∗ (x(tf ), tf ) = 0, where tf is the final time. (25) 1 2 1

E = Ek + Ep ; Ek = Jω ; Ep = Kθ2 (35)
2 2
The optimal cost function J* is quadratic:
where Ek is the kinetic energy and Ep is the potential
energy,
J ∗ (x, t) = xT P (t)x (26)
We can write the cost function J in LQR problem as the
So, total energy of the system.
That means that Q will be:
∂J ∗ (x(t), t)
= Jt∗ = xT Ṗ (t)x  
∂t (27) Jx 0 0 0
∂J ∗ (x(t), t) 1  0 Kx 0 0 
= Jx∗ = 2P (t)x Q=   (36)
∂x 2  0 0 J y 0 
The Hamiltonian: 0 0 0 Ky
For our system, we set:
H = xT Qx + uT Ru + Jx∗T (Ax + Bu)
 
(28) 0.4744 0 0 0
 0 0.05 0 0 
To solve the problem, we have to minimize the Hamiltonian Q = 1013 
 0
,
0 0.2285 0 
with respect to u:
  0 0 0 0.05
1 0
∂H R=
= 2Ru + (Jx∗T B)T = 2Ru + B T Jx∗ = 0 (29) 0 1
∂u and we obtain:
which gives
 
1 0.0059 0.0071 0 −0.0001
u = − R−1 B T Jx∗ (30) K= (37)
0.0001 0.0001 0.0028 0.0071
2
After calculations and replacements, the HJB equation be-
 
0.0071 −0.0001
comes: N= (38)
0.0001 0.0071
The simulink model shown in Fig.6 is used to simulate the
xT Ṗ (t)x + X T Qx + (−R−1 B T P x)T R(−R−1 B T P x) behavior of the closed loop system.

618
 Fig. 6: Simulink Model for feedback control

The step response can be visualized in Fig. 7 where it can

be observed that the settling time is much bigger than the one Fig. 8: Closed-loop step response using modified LQR
obtained with the pole-placement method. weight matrices

V. E XPERIMENTAL RESULTS

Based on the scenarios simulated and presented in the

previous chapter, for experimental test (39) and (40) are used.
As presented in the first chapter, the motors that drive the axes
are controlled with MC206X Motion Coordinator. This have
five servo gains [15]: Proportional gain: P Gain, Integral
gain: I Gain, Derivative gain: D Gain, Output velocity gain:
OV Gain, Velocity Feed-Forward Gain: V F F Gain. The
block diagram that represents the closed loop control system
that can be implemented is represented in Fig. 9.

Fig. 7: Closed-loop step response using LQR

In order to obtain some results that are closer to the one

obtained with the first method, several testes were made and
some conclusions can be mentioned: in order to improve
the response of the closed loop system using the energy
LQR method presented above, it is necessary to change both
R matrix and the weights that correspond to the position
variable. The best results,which are almost the same as the one
obtained with pole placement method, were obtained using the
following
 matrices: 
0.0005 0 0 0
 0 0.85 0 0 
Q= , Fig. 9: Servo loop diagram
 0 0 0.0002 0 
0 0 0 0.25
 
0.01 0 For implementation, we take into consideration the first two
R=
0 0.01 elements from matrix K for X-axis and last two elements for Y-
axis, because those are the one that correspond to the velocity
The simulated step response using the following control and position gain for each axis.
matrices:
So, we set: For X-axis: P Gain = 9.21, OV Gain = 0.21,
and for Y-axis: P Gain = 4.55, OV Gain = 0.13;

0.0059 0.0071 0 −0.0001
 During the experimental tests, several cases were taken into
K= (39) consideration: following a circle path with different radius
0.0001 0.0001 0.0028 0.0071
or squares paths with different size. The resultant path was
reproduced based on the signals and was represented, in each
case, compared to the ones that were obtained using PG ains
 
0.0071 −0.0001
N= (40) set after several tries. [16]
0.0001 0.0071

619
 reason for the error occurrence is the ball screw translational
axis construction and the friction that appears.
In this paper, two methods for designing a controller for an
small size CNC system are presented. Both of the methods
are good for a MIMO system and both of them are easy to
implement. For the LQR problem, an energy based algorithm
of choosing the weight matrices was proposed. At the first
simulation, the controller seemed to work good but slower
than the one obtained with pole placement algorithm. Starting
from the parameters obtained using energy based matrix, the
weights were modified in order to obtain almost the same
result as in the first method. So, the energy based LQR seems
to be a good technique in terms of control and stability, but if
more precise performances is needed, several tests are needed.
So, based on this, we can emphasis that choosing Q and R is
Fig. 10: Results with and without using the designed part art, part science.
controller-squares drawn As a future development, other tests and control algorithms
have to be tried in order to obtain a better fitting. Also,
friction models and compensation should be done for the same
purpose.
R EFERENCES
[1] Kücük Hüseyin Koc, Emine Seda Erdinler, Ender Hazir, and Emel
Öztürk. Effect of cnc application parameters on wooden surface quality.
Measurement, 107:12–18, 2017.
[2] Armin Afkhamifar, Dario Antonelli, and Paolo Chiabert. Variational
analysis for cnc milling process. Procedia CIRP, 43:118–123, 2016.
[3] Xianfeng Chen and Weiming Zhang. Research on the interpolation
technology of cnc system. Procedia engineering, 174:1252–1256, 2017.
[4] Zhiqian Sang and Xun Xu. The framework of a cloud-based cnc system.
Procedia CIRP, 63:82–88, 2017.
[5] SK Jha et al. Comparative study of different classical and modern control
techniques for the position control of sophisticated mechatronic system.
Procedia Computer Science, 93:1038–1045, 2016.
[6] Zhen-yuan Jia, Jian-wei Ma, De-ning Song, Fu-ji Wang, and Wei
Liu. A review of contouring-error reduction method in multi-axis cnc
machining. International Journal of Machine Tools and Manufacture,
2017.
Fig. 11: Results with and without using the designed [7] Kaan Erkorkmaz and Wilson Wong. Rapid identification technique
controller-circles drawn for virtual cnc drives. International Journal of Machine Tools and
Manufacture, 47(9):1381–1392, 2007.
[8] Yoram Koren. Cross-coupled biaxial computer control for manufacturing
In the figures above, it can be seen that in each case a small systems. Journal of Dynamic Systems, Measurement, and Control,
error exists. For each case, a fit performance parameter was 102(4):265–272, 1980.
[9] Dumitru Popescu, Amira Gharbi, Dan Stefanoiu, and Pierre Borne.
calculated and the result were represented in Table.1. Linear identification of closed-loop systems. Process Control Design
for Industrial Applications, pages 21–92, 2017.
Using designed controller No designed controller [10] Lennart Ljung. System identification. In Signal analysis and prediction,
X Axis Y Axis X Axis Y Axis pages 163–173. Springer, 1998.
Circle R=1cm 93.2% 96% 92% 92.7% [11] Honghai Wang, Jianchang Liu, Xia Yu, Shubin Tan, and Yu Zhang.
R=2cm 95% 98% 93% 97% New results on pid controller design of discrete-time systems via pole
Square L=3cm 80% 80% 79% 76% placement. IFAC-PapersOnLine, 50(1):6703–6708, 2017.
L=5cm 84.5% 86% 82% 84.1% [12] Stefano Longo, Eric C Kerrigan, and George A Constantinides.
Constrained lqr for low-precision data representation. Automatica,
50(1):162–168, 2014.
TABLE I: Fit parameter calculated for each case [13] Loı̈c Bourdin and Emmanuel Trélat. Linear–quadratic optimal sampled-
data control problems: Convergence result and riccati theory. Automat-
ica, 79:273–281, 2017.
VI. C ONCLUSION [14] João Marcos Kanieski, Emerson Giovani Carati, and Rafael Cardoso.
An energy based lqr tuning approach applied for uninterruptible power
As a conclusion, based on the latest results presented in supplies. In Circuits and Systems (LASCAS), 2010 First IEEE Latin
the end of the previous chapter, we can say that the designed American Symposium on, pages 41–44. IEEE, 2010.
controller gives a better precision in terms of following a path. [15] Trio motion technology motion coordinator technical reference manual.
2008.
However, we also observed that the error is bigger in case [16] Mariana Ratiu and Alexandru Rus. Modeling of the trajectory-generating
of square path and that could be because at the corners of equipments. In Engineering of Modern Electric Systems (EMES), 2017
the square the acceleration and deceleration are big. Also, a 14th International Conference on, pages 216–219. IEEE, 2017.

620