Review of DC Motor Modeling and Linear Control: Theory with Laboratory Tests

Abstract

This review paper introduces the modeling, measurement, identification and control of direct current motors based on the state space modeling and the transfer function representation. These models are identified by real laboratory measurements, and the simulated results are compared with the measurements. Continuous-time and discrete-time PID (Proportional-Integral-Derivative) controllers, discrete-time state feedback and linear quadratic controllers are designed mathematically. The designed controllers are realized by the microcontroller Arduino UNO, and the behavior of the controllers is compared and analyzed. The noisy current signal has been measured by a discrete-time observer, steady-state Kalman filtering is also studied. The practical results of the implemented controllers support the theoretical results very well.

1. Introduction

Direct current motors (DC motors) are a well-known type of electric motor that convert electrical energy into mechanical energy. DC motors find many applications across a wide range of industries and fields due to their ease of use, simplicity, controllability, and versatility. These motors can be found in various automotive applications, industrial conveyor belt systems, robotic arms, CNC machines, and various actuating and power systems, as well as in medical applications or in everyday devices like home appliances, consumer electronics and so forth. The literature on electrical machinery is vast, see, e.g., [1,2,3,4,5,6].

DC motors operate on the physical principle of Lorentz force, where a current-carrying conductor placed in a magnetic field experiences a force: $\overrightarrow{\mathbf{F}}=I\overrightarrow{\mathbf{dl}}\times \overrightarrow{\mathbf{B}}$, where I is the DC current flowing through the conductor part with the infinitesimal length of $\overrightarrow{\mathbf{dl}}$, and $\overrightarrow{\mathbf{B}}$ is the magnetic flux density [7,8]. This force results in the torque that causes the rotating part (rotor) to turn, resulting in mechanical motion.

The homogeneous magnetic field can be generated by windings supplied by DC current or by magnets arranged around the rotating part. The part of the motor where these windings or the magnets are placed is called the stator [1,2,3].

A frame made of wire is placed in a homogeneous magnetic field as shown in Figure 1 [9]. The current I is flowing through the wire of the frame. In this special case, when the vectors $\overrightarrow{\mathbf{dl}}$ and $\overrightarrow{\mathbf{B}}$ are perpendicular to each other, the Lorentz force equation has the simple form of $F=BIl$, where l is the side length of the wireframe, and $B=|\overrightarrow{\mathbf{B}}|$ is the constant magnetic flux density amplitude. The Lorentz force is constant, and perpendicular to the magnetic field and the frame side, too.

Then, a torque T is acting on the frame: $T=Fdcos\vartheta =BIldcos\vartheta$, where d is the distance between the two sides of the frame, and $\vartheta =\vartheta \left(t\right)$ is the angle between the frame and the horizontal ($\vartheta =0$ in Figure 1). Since the frame can rotate, meaning that the angle can change, consequently the angle is a function of time.

Figure 1 shows only one frame for simplicity. The rotor in a real DC motor consists of a number of frames, which are rotated relative to each other and all frames are connected to the so-called commutator [1,2,3,9].

Since the motor has several frames [1,2,3,9], thanks to the commutators, $cos\vartheta \approx 1$, because there is usually a frame with approximately the maximum value of the torque, and $T=KI$ can be used as an approximation. In other words, torque is proportional to the current, and K is the motor constant depending on the magnetic flux and the motor dimensions (here $K=Bld$).

The frame is rotated by the above-mentioned torque. The voltage $u_{\mathrm{i}}=u_{\mathrm{i}}\left(t\right)$ is induced according to Lenz law, i.e., to the motion: $u_{\mathrm{i}}=-\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}t}$, where $\mathrm{\Phi }=Bldsin\vartheta$ is the flux of one frame, i.e., $u_{\mathrm{i}}=-Bld\dot{\vartheta }cos\vartheta$. Since the motor is built up by several frames, $u_{\mathrm{i}}=-K\omega$ can be supposed, and $\omega =\dot{\vartheta }$ is the angular velocity of the rotor.

The motor constant K is supposed to be a constant value in this paper. It is a motor parameter identified by measurements.

As mentioned, the DC motor consists of several rotating frames [1,2,3,9], and only one frame acts at a moment. The switching between the frames can be realized by the commutator.

There are different types of DC motors [1,2,3,10], including brushed DC motors and brushless DC motors (BLDC motors). Brushed motors have brushes and a commutator to control the direction of current flow in the motor windings. Simple and economical brushed DC motors require maintenance due to brush wear and generate more electrical noise. Brushless motors use electronic commutation instead of brushes. They typically have a higher efficiency, longer lifespan, and lower maintenance requirements compared to brushed DC motors.

DC motors can be controlled in various ways, including voltage control, pulse-width modulation (PWM), or by using specialized motor controllers or drivers. These control methods allow for precise speed and torque control and precise positioning. Generally, feedback control systems use sensors to measure parameters such as motor speed, torque, or position and provide feedback to the motor controller to adjust the operation accordingly in a closed-loop system. Feedback control enables the precise control of the motor performance and allows it to maintain the desired speed, torque, or position even in the presence of disturbances or changes in operating conditions. The literature on drive is also vast, e.g., [4,5,6,10].

To design control systems [11,12,13,14], model based approach is the modern way [15,16,17,18,19,20,21], which is followed in this paper. The studied controllers are the PID controller family [9,17,18,19,20,21,22,23,24,25,26,27,28,29] and the state feedback controllers [10,16,18,19,21,30,31]. PID controller is a type of feedback control system widely used in industrial and automation applications to regulate the behavior of dynamic systems. PID stands for proportional, integral, and derivative, which are the three terms that make up the control action. The controller gains can be set up by the phase margin, for example. State feedback control is a modern technique used in control. The behavior of a dynamic system is controlled by applying feedback based on the state variables of the system, i.e., state feedback control uses measurements of the system’s internal state variables. Various control design techniques, such as pole placement, linear quadratic regulator, and optimal control methods can be used to determine the appropriate feedback gains. The so-called state observer is also studied in this work, which is used to estimate the internal state variables of the dynamic system based on available and deterministic input and output measurements, i.e., an observer estimates the state variables based on the system inputs and outputs. The Kalman filter as an optimal state estimator is also implemented to estimate the state of a dynamic system based on noisy measurements of the system’s stochastic input and output. Here, only linear control methods and full-order observers/estimators are shown.

Modeling the dynamic behavior of a DC motor involves describing its electrical and mechanical characteristics mathematically [9,17,20,21,31,32,33,34,35,36,37,38,39,40]. The dynamic model typically considers the interaction between electrical quantities (voltage and current) and mechanical quantities (torque, speed, and position) in the form of a system of coupled differential equations in the time domain or applying the transfer function between the input and output of the system in the frequency domain or in the Laplace transform domain. Differential equation representation is also called the state space model. Linear and nonlinear models are also shown in the literature, here only linear modeling is studied, and all the nonlinearities are neglected. Major nonlinearity is the dead-zone type friction [29,35,36,37,39,41,42,43].

The dynamic model will be described step by step in Section 2. State feedback controller design is based on the state space representation of the dynamic system, while PID controllers are designed by using the transfer function. PID controllers and state feedback-based controllers are presented and designed separately in this work.

After motor modeling, the system of differential equations or the transfer function can be set up by system identification. The offline (batch method) or online (real-time method) least square parameter estimation are core methods to model identification [30,32,33,35,44,45,46]. The used measurement data and the steps of identification are detailed in Section 2, where the motor parameters (armature resistance and inductance, motor constant, moment of inertia, viscous friction constant) are identified directly from measurements based on the method shown in [20,30,31,34,39,47,48,49]. It will be proved that the mathematical model of the studied DC motor is linear so linear system modeling, identification and linear control methods can be studied in a comprehensive way. It is supposed that the motor parameters are constant and do not vary with time or temperature.

Control design can be performed by using continuous-time model and discrete-time model [11,12,13,14]. First, the continuous-time model and controller design are presented, and then the discrete-time model and design are studied in Section 2. Processing sensor data, calculating control signal, and driving the actuator to achieve desired system behavior are realized by the microcontroller Arduino UNO [50]. Finally, a simple, beginner-friendly and cheap microcontroller environment has been built to study DC motor measurements, modeling, identification and control. The implementation details are shown in Section 3.

The paper focuses on feedback control of the motor speed, torque control can not be performed by the realized test bench.

The built laboratory equipment is suitable for the comprehensive testing of small DC motors, for building and identifying the motor model, as well as for testing the designed controllers. The equipment and the summarized theoretical results can be used for educational purposes too. The aim of the paper is that the measurements with the presented test bench as well as the controller design results can be reproduced. The measurement setup can be used to test new controller algorithms, here only conventional linear techniques are presented in detail. The novelty of this review paper is the presentation, modeling, analysis, and control of a second-order linear and time-invariant system in a comprehensive way. The results based on a literature search can be applied to other second-order systems; moreover, the results can be reproduced.

Section 2 presents the realized test bench, the theoretical background of state space modeling and transfer function representation, PID controllers and state feedback design, equations of the observer as well as the Kalman filter. The results are presented in Section 3, and Section 4 is for discussion. Some important Arduino UNO codes are in the Appendix A, Appendix B and Appendix C.

2. Materials and Methods

2.1. DC Motor Test Bench

A simple laboratory test bench has been implemented to demonstrate the analysis of the DC motor behavior and the designed controllers. The measurement setup shown in Figure 2 can be used to perform both open-loop tests and closed-loop controls.

The type of permanent magnet low-power motor under test is Transtecno EC020.120 [51]. The nominal power is 20 W, and the nominal speed is 2850 RPM (round-per-minute). The motor speed is measured by a disc (with eight black glued strips evenly spaced at ${45}^{\circ }$) rotating between the transmitter and receiver sides of a Sharp GP1A22LC optocoupler [52]. The speed measurement has been implemented in the microcontroller with an interrupt. The type of Hall-element sensor used to measure the current is ACS 712 [53], which can be easily interfaced with the Arduino controller [50]. Here, the well-known Arduino UNO microcontroller has been installed. Details of the power supply and the electronics controlling the motor can be found in the thesis [54].

The open-loop measurement is presented in Figure 3 means that the Arduino microcontroller controls the motor according to a predefined time-dependent voltage signal $u_{\mathrm{a}}=u_{\mathrm{a}}\left(t\right)$, and the measured current $i_{\mathrm{a}}=i_{\mathrm{a}}\left(t\right)$; the speed waveforms $n=n\left(t\right)$ are only recorded for further use. The measured data are fed to the computer via a serial port. This setup is suitable for testing the motor and taking measurements to be recorded for model identification. The index ‘a’ stands for the armature coil, which is fed by the voltage $u_{\mathrm{a}}$ and the current flowing through the coil is $i_{\mathrm{a}}$. In this permanent magnet motor, the armature coil is wound on the rotor. Open-loop measurements are useful for system identification, because motor model parameters are not provided by the manufacturer, i.e., model parameters must be identified via measurements.

In a closed loop arrangement, i.e., control as known from control theory, the speed and in some cases also the current signals are fed back, as shown in Figure 4. In this work, the main control objective is speed control, so this signal is always fed back. The feedback of the current signal is indicated by a dashed line, highlighting that it is not always used (current is necessary only in the case of state feedback control, for simple PID it is not essential). In the case of control, the microcontroller does not just record the measured signals but uses them for control.

Note, that the microcontroller can define the voltage $u_{\mathrm{a}}$ by a PWM value (i.e., a pulse width modulation signal with the appropriate duty cycle), and the power supply and amplifier block results the voltage in the unit of volt.

2.2. DC Motor Modeling

The DC motor equations can be written using the equivalent circuit diagram shown in Figure 5 [3,4,5,6,9,17,20,21]. The excitation of the armature is the terminal voltage $u_{\mathrm{a}}=u_{\mathrm{a}}\left(t\right)$, also called the armature voltage and the armature current in the rotor winding is denoted by $i_{\mathrm{a}}=i_{\mathrm{a}}\left(t\right)$. The electric block represents the armature circuit, which can be modeled as a simple serial-connected RL circuit, and the voltage $u_{\mathrm{i}}=u_{\mathrm{i}}\left(t\right)$ induced by the rotation of the rotor acts against the armature voltage so that it is connected in series with the RL circuit in the model, thus reducing the voltage on the RL-elements. The negative sign shown in the introduction is neglected by this reference direction. The mechanical block of the model describes the rotation of the rotor and the mechanical characteristics of the rotational motion.

The electrical equation of a DC motor based on the Kirchhoff’s voltage law is as follows:

$$R_{\mathrm{a}}{i}_{\mathrm{a}}+L_{\mathrm{a}}{\dot{i}}_{\mathrm{a}}+u_{\mathrm{i}}=u_{\mathrm{a}},$$

(1)

where $R_{\mathrm{a}}$ and $L_{\mathrm{a}}$ are the resistance and inductance of the serial-connected connected RL circuit of the armature. The induced voltage is proportional to the angular velocity of the rotor $\omega =\omega \left(t\right)$, i.e.,

$$u_{\mathrm{i}}=K\omega ,$$

(2)

where K is the constant specific to the machine. The dot above the letter $i_{\mathrm{a}}$ denotes the time derivative, ${\dot{i}}_{\mathrm{a}}=\frac{\mathrm{d}{i}_{\mathrm{a}}}{\mathrm{d}t}$.

The mechanical equation of the machine is as follows:

$$J\dot{\omega }=T-T_{\mathrm{L}}-k\omega ,$$

(3)

where

$$T=K{i}_{\mathrm{a}}$$

(4)

is the torque proportional to the armature current ($T=T\left(t\right)$), $T_{\mathrm{L}}=T_{\mathrm{L}}\left(t\right)$ is the unknown load torque, J is the moment of inertia of the rotating mass, and k is the friction coefficient. The load torque and friction reduce the torque of the motor accordingly.

The following system of linear ordinary differential equations describes the state space model:

$$\begin{array}{cc}\hfill {\dot{i}}_{\mathrm{a}}=& -\frac{R_{\mathrm{a}}}{L_{\mathrm{a}}}{i}_{\mathrm{a}}-\frac{K}{L_{\mathrm{a}}}\omega +\frac{1}{L_{\mathrm{a}}}{u}_{\mathrm{a}},\hfill \\ \hfill \dot{\omega }=& \frac{K}{J}{i}_{\mathrm{a}}-\frac{k}{J}\omega -\frac{1}{J}{T}_{\mathrm{L}},\hfill \end{array}$$

(5)

or with matrix notation:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{i}_{\mathrm{a}}\\ \omega \end{array}\right]=& \left[\begin{array}{cc}-\frac{R_{\mathrm{a}}}{L_{\mathrm{a}}}& -\frac{K}{L_{\mathrm{a}}}\\ \frac{K}{J}& -\frac{k}{J}\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{a}}\\ \omega \end{array}\right]+\left[\begin{array}{c}\frac{1}{L_{\mathrm{a}}}\\ 0\end{array}\right]u_{\mathrm{a}}+\left[\begin{array}{c}0\\ -\frac{1}{J}\end{array}\right]T_{\mathrm{L}},\hfill \\ \hfill y=& \left[\begin{array}{cc}0& 1\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{a}}\\ \omega \end{array}\right]+0u_{\mathrm{a}}.\hfill \end{array}$$

(6)

The unknown load torque $T_{\mathrm{L}}$ can be considered as a disturbance signal, so the input of the system is the armature voltage $u_{\mathrm{a}}$ and the output is the angular velocity $\omega$. The state variables are the armature current and the angular velocity:

$$\mathbf{x}={\left[\begin{array}{cc}{i}_{\mathrm{a}}& \omega \end{array}\right]}^{\mathrm{T}},$$

(7)

moreover

$$\mathbf{A}=\left[\begin{array}{cc}-\frac{R_{\mathrm{a}}}{L_{\mathrm{a}}}& -\frac{K}{L_{\mathrm{a}}}\\ \frac{K}{J}& -\frac{k}{J}\end{array}\right],\mathbf{b}=\left[\begin{array}{c}\frac{1}{L_{\mathrm{a}}}\\ 0\end{array}\right],{\mathbf{c}}^{\mathrm{T}}=\left[\begin{array}{cc}0& 1\end{array}\right],D=0$$

(8)

are the elements of the general description of the linear state space modeling:

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}=& \mathbf{A}\mathbf{x}+\mathbf{b}u,\hfill \\ \hfill y=& {\mathbf{c}}^{\mathrm{T}}\mathbf{x}+Du.\hfill \end{array}$$

(9)

Note, that the angular velocity $\omega$ is often used instead of the rotation speed n. The unit of angular velocity is $\mathrm{rad}/\mathrm{s}$ and the unit of rotation speed is round-per-minute, abbreviated to RPM, i.e., the number of revolutions per minute, or the number of times the motor shaft turns round in one minute. The relationship between the two is as follows:

$$n\left[\mathrm{RPM}\right]=60\frac{\omega [\mathrm{rad}/\mathrm{s}]}{2\pi },\mathrm{and}\omega [\mathrm{rad}/\mathrm{s}]=2\pi \frac{n\left[\mathrm{RPM}\right]}{60}.$$

(10)

In the physical model $\omega$ is used, n is used in microcontroller programs and in the presentation of the measurement results.

The state space modeling will be used in the state feedback controller design. PID controller design is based on the transfer function of the system, which can be derived as follows.

Laplace transforming the differential equation system (6), the following algebraic equations can be obtained:

$$\begin{array}{cc}\hfill (R_{\mathrm{a}}+s{L}_{\mathrm{a}})I_{\mathrm{a}}=& U_{\mathrm{a}}-K\mathrm{\Omega },\hfill \\ \hfill (k+sJ)\mathrm{\Omega }=& K{I}_{\mathrm{a}}-\mathcal{L}\left\{T_{\mathrm{L}}\right\}.\hfill \end{array}$$

(11)

Here, $I_{\mathrm{a}}=I_{\mathrm{a}}\left(s\right)=\mathcal{L}\left\{i_{\mathrm{a}}\left(t\right)\right\}$, $\mathrm{\Omega }=\mathrm{\Omega }\left(s\right)=\mathcal{L}\left\{\omega \right(t\left)\right\}$, and $U_{\mathrm{a}}=U_{\mathrm{a}}\left(s\right)=\mathcal{L}\left\{u_{\mathrm{a}}\left(t\right)\right\}$ denote the Laplace transform of the armature current, angular velocity and armature voltage. The transformed variable $\mathcal{L}\left\{T_{\mathrm{L}}\right\}$ stands for the Laplace transform of the load torque, used in the equations but it will be omitted at the end.

The Laplace transform can be easily used to draw the block diagram shown in Figure 6. Here, it can be clearly seen that the input of the electric block is the voltage $U_{\mathrm{a}}-K\mathrm{\Omega }$ and the output is the current $I_{\mathrm{a}}$, from the first equation of (11)

$$I_{\mathrm{a}}=\frac{U_{\mathrm{a}}-K\mathrm{\Omega }}{R_{\mathrm{a}}+s{L}_{\mathrm{a}}}.$$

(12)

This current is the input of the mechanical block, whose output is the angular velocity $\mathrm{\Omega }$ and the voltage $U_{\mathrm{i}}$, according to the second equation of (11):

$$\mathrm{\Omega }=\frac{K{I}_{\mathrm{a}}-\mathcal{L}\left\{M_{\mathrm{T}}\right\}}{k+sJ},U_{\mathrm{i}}=K\mathrm{\Omega }.$$

(13)

In the following, the load torque is omitted.

From (13)

$$I_{\mathrm{a}}=\frac{k+sJ}{K}\mathrm{\Omega },$$

(14)

which is substituted into (12) to obtain the transfer function for the relationship between the angular velocity and the armature voltage:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{K}{s^{2}J{L}_{\mathrm{a}}+s(k{L}_{\mathrm{a}}+J{R}_{\mathrm{a}})+k{R}_{\mathrm{a}}+K^2},$$

(15)

so that a transfer function describing the relationship between the armature current and the armature voltage can be written as follows:

$$W_{I_{\mathrm{a}}}=\frac{I_{\mathrm{a}}}{U_{\mathrm{a}}}=\frac{k+sJ}{s^{2}J{L}_{\mathrm{a}}+s(k{L}_{\mathrm{a}}+J{R}_{\mathrm{a}})+k{R}_{\mathrm{a}}+K^2}.$$

(16)

By neglecting the friction coefficient, the model can be simplified, i.e., if $k\approx 0$, then the transfer function of the system is as follows:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{K}{s^{2}J{L}_{\mathrm{a}}+sJ{R}_{\mathrm{a}}+K^2},$$

(17)

and

$$W_{I_{\mathrm{a}}}=\frac{I_{\mathrm{a}}}{U_{\mathrm{a}}}=\frac{sJ}{s^{2}J{L}_{\mathrm{a}}+sJ{R}_{\mathrm{a}}+K^2}.$$

(18)

The former transfer function $W_{\mathrm{\Omega }}$, divided by $K^2$ in both the numerator and the denominator, can be expressed as follows:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{1/K}{s^2\frac{J{R}_{\mathrm{a}}}{K^2}\frac{L_{\mathrm{a}}}{R_{\mathrm{a}}}+s\frac{J{R}_{\mathrm{a}}}{K^2}+1}=\frac{1/K}{s^2{\tau }_{\mathrm{m}}{\tau }_{\mathrm{e}}+s{\tau }_{\mathrm{m}}+1},$$

(19)

introduced by

$${\tau }_{\mathrm{e}}=\frac{L_{\mathrm{a}}}{R_{\mathrm{a}}},\mathrm{and}{\tau }_{\mathrm{m}}=\frac{J{R}_{\mathrm{a}}}{K^2},$$

(20)

the so-called electric and mechanical time constants.

For identification and control design purposes, the general transfer function of the second-order system can be written independently of the physical background:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{\alpha }{s^2+\delta s+\epsilon },\mathrm{and}{W}_{I_{\mathrm{a}}}=\frac{I_{\mathrm{a}}}{U_{\mathrm{a}}}=\frac{\beta s+\gamma }{s^2+\delta s+\epsilon },$$

(21)

where $\alpha$, $\beta$, $\gamma$, $\delta$ and $\epsilon$ are parameters.

Identification is detailed in Section 3.2.

2.3. Linear Controller Design

Controller design can be performed after model selection and model parameter identification.

The classical PID controllers can be designed by using the transfer function which output is the motor speed. Here, the design steps based on the phase margin will be followed. The modern state feedback controllers are based on the pole placement method, which can be designed by using the Ackermann equation or the framework of linear quadratic regulation.

These methods are presented next.

2.3.1. PID Controllers

To design the PID controllers, the block diagram shown in Figure 7 is used [11,12,13,18,21]. Here, $W_{\mathrm{P}}=W_{\mathrm{\Omega }}$ denotes the transfer function of the plant to be controlled and $W_{\mathrm{C}}$ denotes the transfer function of the controller to be designed.

The controller input $\Delta \omega$ is the difference between the reference signal ${\omega }_{\mathrm{ref}}$ and the output signal $\omega$ of the plant: $\Delta \omega ={\omega }_{\mathrm{ref}}-\omega$, and the output is the control signal $u_{\mathrm{a}}$, which is the armature voltage driving the motor.

The signals presented in Figure 7 can be generated using the transfer functions. Let the open loop transfer function

$$W_{\mathrm{O}}=W_{\mathrm{C}}{W}_{\mathrm{P}},$$

(22)

then

$$\Delta \mathrm{\Omega }=\frac{1}{1+W_{\mathrm{O}}}{\mathrm{\Omega }}_{\mathrm{ref}},U_{\mathrm{a}}=\frac{W_{\mathrm{C}}}{1+W_{\mathrm{O}}}{\mathrm{\Omega }}_{\mathrm{ref}},\mathrm{\Omega }=\frac{W_{\mathrm{O}}}{1+W_{\mathrm{O}}}{\mathrm{\Omega }}_{\mathrm{ref}}$$

(23)

are the Laplace transform of the error signal, the control signal, and the output signal, and ${\mathrm{\Omega }}_{\mathrm{ref}}$ is the Laplace transform of the reference signal.

The transfer function of the plant must be transformed into a multiplicative form using the time constant, i.e., from (21)

$$W_{\mathrm{P}}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{\alpha }{s^2+\delta s+\epsilon }=\frac{\alpha }{(s+p_1)(s+p_2)}=\frac{\alpha /p_1{p}_2}{(1+s{\tau }_1)(1+s{\tau }_2)},$$

(24)

where ${\tau }_1$ and ${\tau }_2$ are the time constants of the plant. The motor usually has two stable real poles, $p_1<0$ and $p_2<0$, with which ${\tau }_1=-1/p_1$ and ${\tau }_2=-1/p_2$, furthermore, ${\tau }_1>{\tau }_2$.

After this, the continuous-time PID controllers are defined as follows:

• The transfer function of the P controller is constant:

  $$W_{\mathrm{C}}=K_{\mathrm{P}}.$$

  (25)
• The transfer function of the PI controller contains a gain $K_{\mathrm{PI}}$ and a time constant $T_{\mathrm{I}}$ in the following formula:

  $$W_{\mathrm{C}}=K_{\mathrm{PI}}\frac{1+s{T}_{\mathrm{I}}}{s{T}_{\mathrm{I}}},\mathrm{and}{T}_{\mathrm{I}}={\tau }_1.$$

  (26)
• The transfer function of the PD controller is with a gain $K_{\mathrm{PD}}$ and two time constants $T_{\mathrm{D}}$ and $T_{\mathrm{D}}^{\prime }$:

  $$W_{\mathrm{C}}=K_{\mathrm{PD}}\frac{1+s{T}_{\mathrm{D}}}{1+s{T}_{\mathrm{D}}^{\prime }},\mathrm{and}{T}_{\mathrm{D}}={\tau }_2,T_{\mathrm{D}}^{\prime }={\tau }_2/n,$$

  (27)

  where $n>1$ is a controller parameter.
• The transfer function of the applied PID controller is as follows:

  $$W_{\mathrm{C}}=K_{\mathrm{PID}}\frac{1+s{T}_{\mathrm{I}}}{s{T}_{\mathrm{I}}}\frac{1+s{T}_{\mathrm{D}}}{1+s{T}_{\mathrm{D}}^{\prime }},\mathrm{and}{T}_{\mathrm{I}}={\tau }_1,T_{\mathrm{D}}={\tau }_2,T_{\mathrm{D}}^{\prime }={\tau }_2/n.$$

  (28)
  This latter controller is also called a PI-PD controller.

The discrete-time PID controllers are defined as follows:

• The P controller is a gain:

  $$W_{\mathrm{C}}=K_{\mathrm{P}}.$$

  (29)
• The transfer function of the PI controller consists of a gain $K_{\mathrm{PI}}$ and a zero $z_{\mathrm{I}}$ defined by the formula:

  $$W_{\mathrm{C}}=K_{\mathrm{PI}}\frac{z-z_{\mathrm{I}}}{z-1},\mathrm{where}{z}_{\mathrm{I}}={\mathrm{e}}^{-T_{\mathrm{s}}/{\tau }_1}.$$

  (30)
• The transfer function of the PD controller contains a gain $K_{\mathrm{PD}}$ and a zero $z_{\mathrm{D}}$ defined by:

  $$W_{\mathrm{C}}=K_{\mathrm{PD}}\frac{z-z_{\mathrm{D}}}{z},\mathrm{where}{z}_{\mathrm{D}}={\mathrm{e}}^{-T_{\mathrm{s}}/{\tau }_2}.$$

  (31)
• The transfer function of the applied PID controller is as follows:

  $$W_{\mathrm{C}}=K_{\mathrm{PID}}\frac{z-z_{\mathrm{I}}}{z-1}\frac{z-z_{\mathrm{D}}}{z},\mathrm{where}{z}_{\mathrm{I}}={\mathrm{e}}^{-T_{\mathrm{s}}/{\tau }_1},z_{\mathrm{D}}={\mathrm{e}}^{-T_{\mathrm{s}}/{\tau }_2}.$$

  (32)

The time constants ${\tau }_1$ and ${\tau }_2$ are from the continuous-time model of the motor, used to define the zeros $z_{\mathrm{I}}$ and $z_{\mathrm{D}}$ of the discrete-time controllers. Note, however, that the controller gains of the continuous-time model and the discrete-time model are denoted by the same notation, of course, their values are not the same.

The sampling time period is denoted by $T_{\mathrm{s}}$. The sampling frequency $f_{\mathrm{s}}=1/T_{\mathrm{s}}$ can be defined in the microcontroller program.

The details on how to design the controller parameters with the help of the phase margin method are shown in Section 3.3.1.

2.3.2. Observer Design and Kalman Filter

In the following, we turn to the design of state-space-based controllers [11,12,13,18,21]. These methods make use of the state vector $\mathbf{x}={\left[\begin{array}{cc}{i}_{\mathrm{a}}& \omega \end{array}\right]}^{\mathrm{T}}$, i.e., the armature current $i_{\mathrm{a}}$ and the angular velocity $\omega$. The latter is also the output of the system.

The signal from the current sensor is rather noisy; however, it will be used only to verify the modeling. In many practical situations, the current signal is not available from measurements.

If a state variable cannot be measured or is deliberately not intended to be measured, the observer can be of good service in determining the deterministic signal of the state variable. In the stochastic case with noise, the state estimate given by the Kalman filter can also be used. In the following, both solutions are considered, and both are presented shortly.

To design the observer or the estimator and the state feedback controller, it is necessary to discretize the continuous-time model in the state variable description (6). The step response equivalent sampling model can be generated using the Matlab R2015a function $\mathtt{c}\mathtt{2}\mathtt{d}$:

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}=& \mathbf{A}\mathbf{x}+\mathbf{b}u,\hfill \\ \hfill y=& {\mathbf{c}}^{\mathrm{T}}\mathbf{x},\hfill \end{array}\to \begin{array}{cc}\hfill \mathbf{x}[k+1]=& 𝚽\mathbf{x}\left[k\right]+𝚪u\left[k\right],\hfill \\ \hfill y\left[k\right]=& {\mathbf{c}}^{\mathrm{T}}\mathbf{x}\left[k\right],\hfill \end{array}$$

(33)

where $𝚽$ is the system matrix, $𝚪$ and ${\mathbf{c}}^{\mathrm{T}}$ are the input and output vectors of the discrete-time model, respectively. The latter is equal to the output vector of the continuous-time model. The system matrix and the input vector depend on the sampling period time $T_{\mathrm{s}}$. The model generates the $(k+1)$-th value of the state vector from the k-th value of the state vector and the input signal, and the k-th value of the system output signal from the k-th value of the state vector.

The differential equation for a continuous-time observer is as follows:

$$\frac{\mathrm{d}}{\mathrm{d}t}\widehat{\mathbf{x}}=\mathbf{F}\widehat{\mathbf{x}}+\mathbf{G}y+\mathbf{H}u,$$

(34)

where

$$\mathbf{F}=\mathbf{A}-\mathbf{G}{\mathbf{c}}^{\mathrm{T}},\mathbf{H}=\mathbf{b},$$

(35)

after the vector $\mathbf{G}$ is designed with the pole placement method to obtain a fast and stable observer.

A continuous-time observer is not used in the following, because solving the differential Equation (34) on a microcontroller is rather time-consuming, instead the discrete-time observer has been implemented. The continuous-time observer will only be used to validate the results of the discrete-time observer.

A continuous-time observer can be designed if the plant is observable, i.e., if the observability matrix

$${\mathbf{M}}_{\mathrm{o}}=\left[\begin{array}{c}{\mathbf{c}}^{\mathrm{T}}\\ {\mathbf{c}}^{\mathrm{T}}\mathbf{A}\\ {\mathbf{c}}^{\mathrm{T}}{\mathbf{A}}^2\\ ⋮\\ {\mathbf{c}}^{\mathrm{T}}{\mathbf{A}}^{n-1}\end{array}\right]$$

(36)

is invertible, and n is the number of the states. Now $n=2$.

In the implementation, the discrete-time observer shown in Figure 8 has been realized. The input of the block representing the DC motor contains the holding unit and the output contains the sampling unit. The observer has two inputs, the deterministic $u\left[k\right]$ and $y\left[k\right]$, and the output is the estimation of the state vector $\widehat{\mathbf{x}}[k+1]$, meaning that the measurement of the states is not necessary, because the state can be computed.

The difference equation for the discrete-time observer is as follows:

$$\widehat{\mathbf{x}}[k+1]=\mathbf{F}\widehat{\mathbf{x}}\left[k\right]+\mathbf{G}y\left[k\right]+\mathbf{H}u\left[k\right],$$

(37)

where

$$\mathbf{F}=𝚽-\mathbf{G}{\mathbf{c}}^{\mathrm{T}},\mathbf{H}=𝚪,$$

(38)

after the vector $\mathbf{G}$ is designed with the pole placement method to obtain a fast and stable observer. Pole placement can be designed by the Ackermann’s formula:

$$\mathbf{G}={\left[\begin{array}{ccccc}0& 0& \cdots & 0& 1\end{array}\right]}^{\mathrm{T}}{\left({\mathbf{M}}_{\mathrm{o}}^{-1}\right)}^{\mathrm{T}}{\phi }_{\mathrm{o}}\left({𝚽}^{\mathrm{T}}\right),$$

(39)

where ${\mathbf{M}}_{\mathrm{o}}$ is the observability matrix, defined next, ${\phi }_{\mathrm{o}}\left(\lambda \right)$ is the characteristic polynomial of the observer with the desired observer eigenvalues ${\lambda }_{\mathrm{o},i}$ ($i=1,\cdots ,n$),

$${\phi }_{\mathrm{o}}\left(\lambda \right)=(\lambda -{\lambda }_{\mathrm{o},1})(\lambda -{\lambda }_{\mathrm{o},2})\cdots (\lambda -{\lambda }_{\mathrm{o},n}).$$

(40)

A discrete-time observer can be designed if the discrete plant model is observable, i.e., if the discrete-time observability matrix

$${\mathbf{M}}_{\mathrm{o}}=\left[\begin{array}{c}{\mathbf{c}}^{\mathrm{T}}\\ {\mathbf{c}}^{\mathrm{T}}𝚽\\ {\mathbf{c}}^{\mathrm{T}}{𝚽}^2\\ ⋮\\ {\mathbf{c}}^{\mathrm{T}}{𝚽}^{n-1}\end{array}\right]$$

(41)

is invertible, and n is the number of the states. Now, $n=2$.

The design of the Kalman filter, suppose that the discrete-time model of the plant is given by the following equations instead of (33):

$$\begin{array}{cc}\hfill \mathbf{x}[k+1]=& 𝚽\mathbf{x}\left[k\right]+𝚪u\left[k\right]+\mathbf{v}\left[k\right],\hfill \\ \hfill y\left[k\right]=& {\mathbf{c}}^{\mathrm{T}}\mathbf{x}\left[k\right]+w\left[k\right],\hfill \end{array}$$

(42)

where the vector $\mathbf{v}\left[k\right]$ and the scalar $w\left[k\right]$ represent the system noise and the output noise, respectively, i.e., the noise model is appended to the state space representation of (33). It is supposed that the noises are Gaussian white noise with zero expected value and the covariances

$${\mathbf{R}}_{\mathrm{v}}=\mathbb{E}\left\{\mathrm{v}\left[k\right]\mathrm{v}{\left[k\right]}^{\mathrm{T}}\right\},R_{\mathrm{w}}=\mathbb{E}\left\{w{\left[k\right]}^2\right\}$$

(43)

are obtained from measurements. All the other covariances are supposed to be zero.

Instead of the iterative algorithm for Kalman filter, the steady-state approximation has been implemented in this work, because the DC motor state space representation is a time invariant model. In the more general former method the Kalman gain is given iteratively, in the latter, the steady-state Kalman gain is given as the solution of the so-called Filter Algebraic Riccati Equation (FARE).

The difference equation of the Kalman filter is the same as that of the discrete-time observer (37); moreover, the matrix $\mathbf{F}$ and the vector $\mathbf{H}$ are calculated in the same way, but the vector $\mathbf{G}$, the so-called steady-state Kalman gain is calculated via the discrete-time FARE:

$$𝚽\mathbf{P}{𝚽}^{\mathrm{T}}-𝚽\mathbf{P}{\mathbf{c}}^{\mathrm{T}}{\left(\mathbf{c}\mathbf{P}{\mathbf{c}}^{\mathrm{T}}+R_{\mathrm{w}}\right)}^{-1}\mathbf{c}\mathbf{P}{𝚽}^{\mathrm{T}}+{\mathbf{R}}_{\mathrm{v}}-\mathbf{P}=\mathbf{0},$$

(44)

from which the steady-state estimation error covariance matrix $\mathbf{P}$ can be obtained, and than the steady-state Kalman gain is constant,

$$\mathbf{G}=𝚽\mathbf{P}{\mathbf{c}}^{\mathrm{T}}{\left(\mathbf{c}\mathbf{P}{\mathbf{c}}^{\mathrm{T}}+R_{\mathrm{w}}\right)}^{-1}.$$

(45)

It is highlighted that the Kalman estimator has two inputs, the stochastic $u\left[k\right]$ and $y\left[k\right]$, and the output is the estimation of the state vector $\widehat{\mathbf{x}}[k+1]$, meaning that the measurement of the states is not necessary, because the state can be estimated.

Section 3.3.2 details the numerical results.

2.3.3. Digital State Feedback Controllers

The block diagram of the system implementing the state feedback control is shown in Figure 9 [11,12,13,18,21]. The input of the motor is the held value of the control signal $u_{\mathrm{a}}\left[k\right]$ with the sampling period time $T_{\mathrm{s}}$, and the output is the sampled angular frequency $\omega \left[k\right]$ or the rotor speed $n\left[k\right]$.

The current signal is not measured directly by the current sensor but is generated by the state observer or by the Kalman filter.

In order to make the speed response to a step reference signal change with zero steady-state error, an integrator is built into the feedforward branch of the closed loop system. The integrator is approximated by the left-side rectangular rule. The integrator input here is $e\left[k\right]=\Delta \omega \left[k\right]$ and the output is denoted by $\delta \left[k\right]$:

$$\frac{1}{s}\cong \frac{T_{\mathrm{s}}{z}^{-1}}{1-z^{-1}}=\frac{T_{\mathrm{s}}}{z-1}=\frac{\mathcal{Z}\left\{\delta \right[k\left]\right\}}{\mathcal{Z}\left\{e\right[k\left]\right\}},$$

(46)

where $\mathcal{Z}\{·\}$ is the z-transform of the signal, and then $\delta [k+1]$ can be written in the following recursive form:

$$\delta [k+1]=\delta \left[k\right]+T_{\mathrm{s}}e\left[k\right].$$

(47)

The appended equation of the state space model extended with the integrator can be written as follows:

$$\left[\begin{array}{c}x[k+1]\\ \delta [k+1]\end{array}\right]=\left[\begin{array}{cc}𝚽& \mathbf{0}\\ T_{\mathrm{s}}{\mathbf{c}}^{\mathrm{T}}& 1\end{array}\right]\left[\begin{array}{c}x\left[k\right]\\ \delta \left[k\right]\end{array}\right]+\left[\begin{array}{c}𝚪\\ 0\end{array}\right]u\left[k\right]={𝚽}^{★}{\mathbf{x}}^{★}\left[k\right]+{𝚪}^{★}u\left[k\right],$$

(48)

for which a state feedback controller can be designed using the Ackermann’s formula:

$$\left[\begin{array}{cc}{\mathbf{k}}^{\mathrm{T}}& k_{\mathrm{I}}\end{array}\right]={\left[\begin{array}{ccccc}0& 0& \cdots & 0& 1\end{array}\right]}^{\mathrm{T}}{\left({\mathbf{M}}_{\mathrm{c}}^{★}\right)}^{-1}{\phi }_{\mathrm{cl}}\left({𝚽}^{★}\right),$$

(49)

where ${\mathbf{M}}_{\mathrm{c}}^{★}$ is the controllability matrix of the appended system, ${\phi }_{\mathrm{cl}}\left(\lambda \right)$ is the characteristic polynomial of the desired closed loop system with the eigenvalues ${\tilde{\lambda }}_i$ ($i=1,\cdots ,n+1$),

$${\phi }_{\mathrm{cl}}\left(\lambda \right)=(\lambda -{\tilde{\lambda }}_1)(\lambda -{\tilde{\lambda }}_2)\cdots (\lambda -{\tilde{\lambda }}_{n+1}).$$

(50)

The eigenvalue for the integrator is also designed, this is the last (($n+1$)-th) eigenvalue. Finally, the controllability matrix of the appended system is defined as follows:

$${\mathbf{M}}_{\mathrm{c}}=\left[\begin{array}{cccccc}{𝚪}^{★}& {𝚽}^{★}{𝚪}^{★}& {\left({𝚽}^{★}\right)}^2{𝚪}^{★}& \cdots & & {\left({𝚽}^{★}\right)}^{n-1}{𝚪}^{★}\end{array}\right].$$

(51)

The system is controllable, if the controllability matrix is invertible, and n is the number of the states. Now, $n=2$.

Linear quadratic control (LQ-control) is a kind of optimal control, which minimizes the function defined by

$$J({\mathbf{x}}^{★}\left[k\right],u\left[k\right])=\frac{1}{2}\sum _{k=0}^{\infty }\left[{\left({\mathbf{x}}^{★}\left[k\right]\right)}^{\mathrm{T}}{\mathbf{Q}}^{★}{\mathbf{x}}^{★}\left[k\right]+R{u}^2\left[k\right]\right],$$

(52)

and the appended system of difference equations of the closed loop system is a linear constraint of the optimization problem.

The first term of the functional provides a stable closed loop system, while the second term can influence the control energy.

Linear quadratic control (LQ-control) scheme also results in the state feedback realization, so the same design steps are applied as in the case of the state feedback control. The feedback vector with which the closed-loop system states and the control signal can minimize the functional calculated differently, as the solution to the Discrete Algebraic Riccati Equation (DARE):

$${\left({𝚽}^{★}\right)}^{\mathrm{T}}{\mathbf{X}}^{★}{𝚽}^{★}-{\left({𝚽}^{★}\right)}^{\mathrm{T}}{\mathbf{X}}^{★}{𝚪}^{★}{\left({\left({𝚪}^{★}\right)}^{\mathrm{T}}{\mathbf{X}}^{★}{𝚪}^{★}+R\right)}^{-1}{\left({𝚪}^{★}\right)}^{\mathrm{T}}{\mathbf{X}}^{★}{𝚽}^{★}+{\mathbf{Q}}^{★}-{\mathbf{X}}^{★}=\mathbf{0},$$

(53)

from which the matrix ${\mathbf{X}}^{★}$ can be obtained, and then the feedback gains are as follows:

$$\left[\begin{array}{cc}{\mathbf{k}}^{\mathrm{T}}& k_{\mathrm{I}}\end{array}\right]={\left({\left({𝚪}^{★}\right)}^{\mathrm{T}}{\mathbf{X}}^{★}{𝚪}^{★}+R\right)}^{-1}{\left({𝚪}^{★}\right)}^{\mathrm{T}}{\mathbf{X}}^{★}{𝚽}^{★}.$$

(54)

Numerical results are shown in Section 3.3.2.

3. Results

3.1. Measuring Some Motor Characteristics

With the open loop arrangement (Figure 3), the static, i.e., the steady-state conditions can be easily investigated. Some data are shown in

Table 1. Measured data in the stationary state.

${\mathit{u}}_{\mathrm{a}}\left[\mathrm{PWM}\right]$	50	60	70	80	90	100	110	120
$u_{\mathrm{a}}\left[\mathrm{V}\right]$	0.50	0.93	1.25	1.76	2.21	2.65	3.08	3.47
$i_{\mathrm{a}}\left[\mathrm{A}\right]$	-	-	-	-	0.04	0.05	0.06	0.08
$n\left[\mathrm{RPM}\right]$	90	245	376	520	650	792	930	1062
$\omega [\mathrm{rad}/\mathrm{s}]$	9.43	25.66	39.38	54.45	68.07	82.94	97.39	111.21
$u_{\mathrm{a}}\left[\mathrm{PWM}\right]$	130	140	150	160	170	180	190	200
$u_{\mathrm{a}}\left[\mathrm{V}\right]$	3.95	4.42	4.89	5.30	5.71	6.12	6.55	6.97
$i_{\mathrm{a}}\left[\mathrm{A}\right]$	0.11	0.12	0.16	0.17	0.19	0.23	0.25	0.26
$n\left[\mathrm{RPM}\right]$	1215	1356	1511	1648	1772	1923	2058	2192
$\omega [\mathrm{rad}/\mathrm{s}]$	127.23	142.00	158.23	172.58	185.56	201.38	215.51	229.55
$u_{\mathrm{a}}\left[\mathrm{PWM}\right]$	210	220	230	240	250	255
$u_{\mathrm{a}}\left[\mathrm{V}\right]$	7.46	7.91	8.31	8.80	9.82	10.4
$i_{\mathrm{a}}\left[\mathrm{A}\right]$	0.28	0.30	0.32	0.35	0.39	0.40
$n\left[\mathrm{RPM}\right]$	2346	2489	2625	2786	3120	3308
$\omega [\mathrm{rad}/\mathrm{s}]$	245.67	260.65	274.89	291.75	326.73	346.41

.

The Arduino UNO uses a PWM signal to adjust the input voltage. The PWM output is 8-bit, i.e., the PWM value can be an integer between 0 and 255, thus determining the duty cycle, i.e., the armature voltage driving the motor.

By exciting the motor with different PWM values ($u_{\mathrm{a}}\left[\mathrm{PWM}\right]$, the first row of Table 1), the steady-state input voltage $u_{\mathrm{a}}\left[\mathrm{V}\right]$, the steady-state armature current $i_{\mathrm{a}}\left[\mathrm{A}\right]$, the rotor speed $n\left[\mathrm{RPM}\right]$ and the angular velocity $\omega [\mathrm{rad}/\mathrm{s}]$ can be measured. The steady state is typically reached within 1 second.

First, the relationship between the integer PWM input value and the stationary motor state variables is examined. In this case, the amplifier circuit is understood to be in the same block as the motor, as highlighted by the dotted line in the block diagram in Figure 3.

Second, the relationship between the input voltage measured in volts and the state variables in the stationary state is discussed. In this second case, the amplifier circuit is understood as a block with the microcontroller, as highlighted by the dashed line in Figure 3.

In the former case, the input-output relationship is nonlinear, while in the latter case, it is linear. In the controller design, the linear approach is used, and this approach is more closely related to the state space description of the motor since its input is the physically interpretable armature voltage measured in volts.

Figure 10 shows one of the static characteristics of the motor, the relationship between the 8-bit PWM excitation signal and the rotor speed. It can be seen that this is nonlinear, and given by the following characteristic:

$$n\left[\mathrm{RPM}\right]=\left\{\begin{array}{ccc}0& \mathrm{if}& u_{\mathrm{a}}\left[\mathrm{PWM}\right]<47,\\ 14.7u_{\mathrm{a}}\left[\mathrm{PWM}\right]-683.3& \mathrm{if}& u_{\mathrm{a}}\left[\mathrm{PWM}\right]\ge 47.\end{array}\right.$$

(55)

This means that the motor shaft will not rotate until the 8-bit PWM excitation reaches (rounded up to the integer value) 47. The small circle indicates the measurement results and the solid line indicates the approximate values according to (55).

All the interpolations performed by linear parameter estimation technique called the least square method [11,12].

In Figure 10, a possible linear approximation is also plotted with a dashed line, with the equation $n\left[\mathrm{RPM}\right]=12.973u_{\mathrm{a}}\left[\mathrm{PWM}\right]$, the slope of which is calculated from the last measured data. This will not be discussed further.

The same type of armature current-PWM excitation relationship is described by

$$i_{\mathrm{a}}\left[\mathrm{A}\right]=\left\{\begin{array}{ccc}0& \mathrm{if}& u_{\mathrm{a}}\left[\mathrm{PWM}\right]<71,\\ 0.002059u_{\mathrm{a}}\left[\mathrm{PWM}\right]-0.14784& \mathrm{if}& u_{\mathrm{a}}\left[\mathrm{PWM}\right]\ge 71.\end{array}\right.$$

(56)

The PWM value of the excitation is related to the voltage measured on the motor terminals by a volt-meter instrument as follows:

$$u_{\mathrm{a}}\left[\mathrm{V}\right]=\left\{\begin{array}{ccc}0& \mathrm{if}& u_{\mathrm{a}}\left[\mathrm{PWM}\right]<43,\\ 0.045521u_{\mathrm{a}}\left[\mathrm{PWM}\right]-1.9336& \mathrm{if}& u_{\mathrm{a}}\left[\mathrm{PWM}\right]\ge 43.\end{array}\right.$$

(57)

The relationship between the armature voltage of the motor measured in volts and the rotor speed is linear according to the measurements shown in Figure 11. The linear approximation is given by

$$n\left[\mathrm{RPM}\right]=318.077u_{\mathrm{a}}\left[\mathrm{V}\right].$$

(58)

The small circle in Figure 11 indicates the measurement results, solid line indicates the approximate values according to the above equations.

The armature current-armature voltage relationship can also be considered linear to a good approximation (see Figure 12):

$$i_{\mathrm{a}}\left[\mathrm{A}\right]=0.03846u_{\mathrm{a}}\left[\mathrm{V}\right].$$

(59)

The characteristic above $0.2\mathrm{A}$ is considered linear, below that the larger deviation is due to the inaccuracy of the sensor, and linear behavior is assumed for the rest of the paper.

3.2. Identification of Motor Parameters

In the steady-state the derivatives vanish in (1)–(4), and the following equations can be obtained:

$$u_{\mathrm{a}}=R_{\mathrm{a}}{i}_{\mathrm{a}}+K\omega ,\mathrm{and}K{i}_{\mathrm{a}}=k\omega ,$$

(60)

from which K and k can be expressed:

$$K=\frac{u_{\mathrm{a}}-R_{\mathrm{a}}{i}_{\mathrm{a}}}{\omega },\mathrm{and}k=\frac{K{i}_{\mathrm{a}}}{\omega }.$$

(61)

It can be seen that the constants K and k depend on the value of the resistance $R_{\mathrm{a}}$.

The block diagram of Figure 6 can be split into two separate parts as shown in Figure 13. The input of the electric part is the voltage $u_{\mathrm{a}}-K\omega$ and the output is the current signal $i_{\mathrm{a}}$. The current is the input of the mechanical part, the output of which is the angular velocity $\omega$. Thus, if the input-output time functions of interest are available from measurements, the unknown coefficients in the blocks can be determined very accurately after a good choice of armature resistance. The inductance of the armature branch $L_{\mathrm{a}}$ and the moment of inertia of the mechanical part J can be identified independently. The identification is, therefore, carried out in the following order:

$$R_{\mathrm{a}}\to (K,k)\to (L_{\mathrm{a}},J).$$

Using the Matlab R2015a $\mathtt{fminbnd}$ function, the motor parameters can be determined using an iterative procedure based on the golden section method as follows:

(1)

Initialization of $R_{\mathrm{a}}$ on the interval $[0.1,\cdots ,10]\mathrm{\Omega }$ according to the rules of the golden section method;

(2)

K and k are defined by (61) as a function of the resistance $R_{\mathrm{a}}$, averaged over the measured stationary results (only the values between $u_{\mathrm{a}}\left[\mathrm{PWM}\right]=150$ and $u_{\mathrm{a}}\left[\mathrm{PWM}\right]=255$ in Table 1 are used);

(3)

Determine the value of the inductance $L_{\mathrm{a}}$ using the function $\mathtt{fminbnd}$ and the differential equation of the electric block shown in Figure 13;

(4)

Determine the value of the moment of inertia J using the function $\mathtt{fminbnd}$ and the differential equation of the mechanical block shown in Figure 13;

(5)

Modify $R_{\mathrm{a}}$ according to the golden section method, and go back to step 2 until the algorithm converges.

In step 3 and in step 4, the following differential equations must be solved:

$${\dot{i}}_{\mathrm{a}}=-\frac{R_{\mathrm{a}}}{L_{\mathrm{a}}}{i}_{\mathrm{a}}+\frac{1}{L_{\mathrm{a}}}(u_{\mathrm{a}}-K\omega ),$$

(62)

and

$$\dot{\omega }=-\frac{k}{\mathrm{\Theta }}\omega +\frac{1}{\mathrm{\Theta }}\left(K{i}_{\mathrm{a}}\right).$$

(63)

In the former case, the time function of the current signal, in the latter the time function of the angular velocity is obtained. The calculated time functions are compared with the measured results and an error is calculated. To conduct this, the motor is excited with a PWM signal of 250, for example, and the time functions of angular velocity and current in the transient are measured. Given the measured data, the inputs of the electric and mechanical block can be generated separately and the differential equation can be solved numerically using one of the Matlab R2015a $\mathtt{ode}$ functions.

The identification results are as follows:

$$\begin{array}{cc}\hfill R_{\mathrm{a}}=& 6.21\mathrm{\Omega },L_{\mathrm{a}}=164\mathrm{mH},\hfill \\ \hfill K=& 0.0234\frac{\mathrm{V}}{\frac{\mathrm{rad}}{\mathrm{s}}},k=2.63·{10}^{-5}\frac{\mathrm{VA}}{{\left(\frac{\mathrm{rad}}{\mathrm{s}}\right)}^2},J=2.31·{10}^{-5}\mathrm{kg}{\mathrm{m}}^2.\hfill \end{array}$$

(64)

Substituting the parameters obtained from the identification in SI units into the state space description (6), the following model can be obtained:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{i}_{\mathrm{a}}\\ \omega \end{array}\right]=\left[\begin{array}{cc}-37.8087& -0.1425\\ 1012.9& -1.136\end{array}\right]\left[\begin{array}{c}{i}_{\mathrm{a}}\\ \omega \end{array}\right]+\left[\begin{array}{c}6.0843\\ 0\end{array}\right]u_{\mathrm{a}}.$$

(65)

Substituting the above-mentioned physical parameters into the relations (15) and (16), the following transfer functions are obtained:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{6176}{s^2+s39+188},\mathrm{and}{W}_{I_{\mathrm{a}}}=\frac{I_{\mathrm{a}}}{U_{\mathrm{a}}}=\frac{s6.1+6.9}{s^2+s39+188}.$$

(66)

For the above state space description, for comparison, the following transfer functions are calculated by the Matlab R2015a function $\mathtt{ss}\mathtt{2}\mathtt{tf}$:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{6162.9}{s^2+s38.9446+187.3085},W_{I_{\mathrm{a}}}=\frac{I_{\mathrm{a}}}{U_{\mathrm{a}}}=\frac{s6.0843,1+6.9115}{s^2+s38.9446+187.3085}.$$

(67)

Another, simpler solution to identify the continuous-time transfer function of the second-order system is as follows.

The parameters in the transfer functions (21) can be determined using the Matlab R2015a $\mathtt{fminsearch}$ function. The error to be minimized is the difference between the angular velocity and current response calculated using the searched transfer function and the measured time functions.

The simulated angular velocity and armature current are obtained as a solution of the following state space model:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=& \left[\begin{array}{cc}-\delta & -\epsilon \\ 1& 0\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]+\left[\begin{array}{c}1\\ 0\end{array}\right]u_{\mathrm{a}},\hfill \\ \hfill y=& {\mathbf{c}}^{\mathrm{T}}\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]\hfill \end{array}$$

(68)

where $x_1$ and $x_2$ are notional state variables and y is the notional output. The  output vector ${\mathbf{c}}^{\mathrm{T}}=\left[0\alpha \right]$ gives the angular velocity, and ${\mathbf{c}}^{\mathrm{T}}=\left[\beta \gamma \right]$ gives the armature current.

To conduct the identification, the motor is excited with a PWM signal of 250, for example, and the time functions of angular velocity and current in the transient are measured.

The identified parameters are as follows:

$$\begin{array}{c}\hfill \alpha =5434,\beta =5.53,\gamma =3.82,\delta =32.5,\epsilon =165,\end{array}$$

(69)

i.e., the transfer functions are as follows:

$$W_{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{5434}{s^2+s32.5+165},\mathrm{and}{W}_{I_{\mathrm{a}}}=\frac{I_{\mathrm{a}}}{U_{\mathrm{a}}}=\frac{s5.53+3.82}{s^2+s32.5+165}.$$

(70)

Figure 14 shows the results of the switch-on transients for different input voltage values. The measured time functions of motor speed and armature current are compared with the model results. For speed, the results are plotted for 100, 150, 200 and 250 8-bit PWM voltage values, for current only the results for 150 and 250 PWM excitations are shown for better visibility. It is clearly visible that the current signal is quite noisy.

The time functions for the 250 PWM value have been used for the identifications.

3.3. Design of Linear Controllers

3.3.1. PID Controllers

The block diagram of the closed-loop control system for speed control used in this study is shown in Figure 15. The plant to be controlled is the DC motor, whose speed $n=n\left(t\right)$ is controlled, and the reference signal is $n_{\mathrm{ref}}=n_{\mathrm{ref}}\left(t\right)$. The difference between the two is the error signal:

$$\Delta n=n_{\mathrm{ref}}-n.$$

(71)

The controller is designed for the transfer function $W_{\mathrm{P}}=W_{\mathrm{\Omega }}=\mathrm{\Omega }/U_{\mathrm{a}}$, i.e., assuming a plant model with input of armature voltage in the unit of volt and with the output of angular velocity in rad/s. The error term must, therefore, be converted from RPM to rad/s, this is conducted by the $n\to \omega$ block. The input of the controller is, therefore, the angular velocity error $\Delta \omega$ and the output is the voltage $u_{\mathrm{a}}$ in the unit of volt. The output of the microcontroller produces a voltage in units of PWM, and the $\mathrm{V}\to \mathrm{PWM}$ block converts the voltage signal of the control algorithm into a PWM value. To conduct this, the relation (57) is rearranged:

$$u_{\mathrm{a}}\left[\mathrm{PWM}\right]=21.96788u_{\mathrm{a}}\left[\mathrm{V}\right]+42.477,$$

(72)

and the result, of course, is rounded to the nearest integer.

The control algorithm also contains a lower and an upper limit below and towards which the control signal must not go. These PWM values are 50 and 250. The motor does not rotate below 50 PWM (50 is used instead of 43 in the context (57)), and the maximum value would be 255, but it is set to 250 in the microcontroller program. These limits should be taken into account when designing and checking the operation of the controller.

Design of the Continuous-Time PID Controllers

The PID controllers can be designed using the method of designing for phase margin using the open loop system Bode diagram. To conduct this, the polynomial-per-polynomial transfer function of the plant is transformed into a multiplicative form using the time constant, i.e.,

$$W_{\mathrm{P}}=\frac{\mathrm{\Omega }}{U_{\mathrm{a}}}=\frac{5434}{s^2+s32.5+165}=\frac{5434}{(s+26.2)(s+6.3)}=\frac{32.9}{(1+s0.159)(1+s0.038)}.$$

(73)

The slow time constant is ${\tau }_1=0.159\mathrm{s}$, and the fast one is ${\tau }_2=0.038\mathrm{s}$, the system poles are $p_1=-6.3\frac{1}{\mathrm{s}}$ and $p_2=-26.2\frac{1}{\mathrm{s}}$, with which ${\tau }_1=-1/p_1$ and ${\tau }_2=-1/p_2$.

As an example, let us design controllers for phase margin ${\phi }_{\mathrm{m}}={60}^{\circ }$. As an illustration, the design of a P-controller using Figure 16 is used. In a Matlab R2015a environment, this can be easily conducted with the function $\mathtt{bode}$. For this:

(1)

Plot the Bode diagram of the transfer function $W_{\mathrm{O}}=W_{\mathrm{C}}{W}_{\mathrm{P}}$ of the open loop with unit gain, i.e., let $K_{\mathrm{P}}=1$ ($K_{\mathrm{PI}}=1$, $K_{\mathrm{PD}}=1$, and $K_{\mathrm{PID}}=1$, respectively), as plotted in Figure 16;

(2)

find the point $-{180}^{\circ }+{\phi }_{\mathrm{m}}$ (which here means the point $-{120}^{\circ }$) on the phase diagram, and read off the cut-off frequency ${\omega }_{\mathrm{c}}$ (in the example ${\omega }_{\mathrm{c}}$ = 23.8 $\mathrm{rad}/\mathrm{s}$ at $-{117}^{\circ }$, of course, the $\mathtt{bode}$ function can be used to obtain a more precise value, the $3^{\circ }$ error is due to the readout given by the graphical display);

(3)

read off the value of the amplitude characteristic $K_{\mathrm{dB}}$ in decibel units at the cut-off frequency (here $K_{\mathrm{dB}}\left({\omega }_{\mathrm{c}}\right)=15.9\mathrm{dB}$);

(4)

determine the controller gain:

$$K_{\mathrm{controller}}={10}^{-K_{\mathrm{dB}}/20},$$

(74)

i.e., $K_{\mathrm{P}}={10}^{-15.9/20}=0.16$, which is a ratio without unit, since the $\mathtt{bode}$ function represents $K_{\mathrm{dB}}$ in decibels according to the relationship of $K_{\mathrm{dB}}=20\lg K_{\mathrm{controller}}$.

Figure 16. Design of P controller for phase margin ${\phi }_{\mathrm{m}}={60}^{\circ }$.

Figure 16. Design of P controller for phase margin ${\phi }_{\mathrm{m}}={60}^{\circ }$.

The design can be conducted without plotting the Bode diagram directly, but then a nonlinear system of equations must be solved:

$$|W_{\mathrm{O}}\left(\mathrm{j}{\omega }_{\mathrm{c}}\right)|=1,-{180}^{\circ }+{\phi }_{\mathrm{m}}=\mathrm{arc}{W}_{\mathrm{O}}\left(\mathrm{j}{\omega }_{\mathrm{c}}\right).$$

(75)

This system of equations can be solved efficiently by calling the $\mathtt{fsolve}$ function. The two unknowns associated with the two equations are the controller gain $K_{\mathrm{controller}}$ and the cut-off frequency ${\omega }_{\mathrm{c}}$.

The step responses of the designed closed-loop control systems are compared in Figure 17, and some characteristic data are summarized in

Table 2. Some data of closed-loop systems for the phase margin ${\phi }_{\mathrm{m}}={60}^{\circ }$.

Controller	${\mathit{\omega }}_{\mathbf{c}}[\mathbf{rad}/\mathbf{s}]$	K	${\mathit{v}}_{max}$	${\mathit{v}}_{\mathbf{st}}$	$\mathbf{\Delta }\mathit{v}[\%]$	${\mathit{u}}_{\mathbf{t}}$	${\mathit{u}}_{\mathbf{st}}$
P	25.28	0.175	0.99	0.85	17.4	6.74	0.026
PI	15.13	0.084	1.09	1.00	9.48	2.87	0.031
PD ($n=2$)	41.78	0.261	1.03	0.89	14.6	19.2	0.027
PI-PD ($n=2$)	41.78	0.261	1.17	1.00	14.9	17.2	0.031
PI-PD ($n=5$)	88.57	0.517	1.12	1.00	12.4	85.1	0.031

. Here, ${\omega }_{\mathrm{c}}$ is the cut-off frequency in units of rad/s, K is the controller gain, $v_{max}$ is the maximum value of the step response, and $v_{\mathrm{st}}$ is the steady-state value of the step response, $\Delta v[\%]=100\frac{v_{max}-v_{\mathrm{st}}}{v_{\mathrm{st}}}$ is the overshoot, $u_{\mathrm{t}}=\frac{u\left(0\right)}{u_{\mathrm{st}}}$ is the overcontrol, and finally $u_{\mathrm{st}}$ is the steady-state value of the control signal.

Table 2 shows that the higher the value of ${\omega }_{\mathrm{c}}$, the faster the control loop is, as can be seen in Figure 17, and that increasing n further increases the controller speed. Also, the higher the value of the loop gain, the faster the controller. Further increasing K will speed up the control loop. However, the acceleration leads to oscillatory behavior, as can be clearly seen in the simulation results, and in fact, large gains can make the closed-loop system unstable. However, fast control is possible with large control signals, as can be seen from the large overcontrol values $u_{\mathrm{t}}$. A trade-off between the two must be found, i.e., to achieve the fastest control possible, but without saturating the actuating device. From the steady-state value $v_{\mathrm{st}}$ it is clear that a controller including an integrator allows accurate operation for the step reference signal. Overshooting occurs in all cases, which can be reduced by increasing the phase margin, reducing the gain, or reducing n. This of course makes the control slower.

The closed loop system behavior is not very accurate when applying only a P-controller. It can be stated that the plant response can be accelerated by a PD block, the steady-state accuracy can be improved by a PI block, while PID altogether results in a fast and accurate control loop.

Discrete-Time Realization of the Continuous-Time PID Controllers

The sampled realization of the above-mentioned continuous-time PID controllers can be conducted by the step response equivalent transformation or by the unconditionally stable backward differential scheme, the right-hand-side rectangle rule, and the Tustin formula. The sampling period time is denoted by $T_{\mathrm{s}}$, i.e., the frequency of the sampling provided by the microcontroller is $f_{\mathrm{s}}=1/T_{\mathrm{s}}$.

The constant gain of the discrete-time P controller is of course the same as the gain of the continuous-time controller.

By using $s=\frac{z-1}{T_{\mathrm{s}}z}$ in the transfer function of the PI controller,

$$W_{\mathrm{C}}={\left.K_{\mathrm{PI}}\frac{1+s{T}_{\mathrm{I}}}{s{T}_{\mathrm{I}}}\right|}_{s=\frac{z-1}{T_{\mathrm{s}}z}}=K_{\mathrm{PI}}\frac{1+\frac{z-1}{T_{\mathrm{s}}z}{T}_{\mathrm{I}}}{\frac{z-1}{T_{\mathrm{s}}z}{T}_{\mathrm{I}}}$$

(76)

can be written. After a short transformation of $W_{\mathrm{C}}\left(z\right)=U\left(z\right)/E\left(z\right)$, the following difference equation can be obtained to obtain the discrete-time control signal:

$$u\left[k\right]=u[k-1]+K_{\mathrm{PI}}\left(1+\frac{T_{\mathrm{s}}}{T_{\mathrm{I}}}\right)e\left[k\right]-K_{\mathrm{PI}}e[k-1].$$

(77)

Using the Tustin formula $s=\frac{2}{T_{\mathrm{s}}}\frac{z-1}{z+1}$ results in:

$$u\left[k\right]=u[k-1]+K_{\mathrm{PI}}\left(1+\frac{T_{\mathrm{s}}}{2T_{\mathrm{I}}}\right)e\left[k\right]-K_{\mathrm{PI}}\left(1-\frac{T_{\mathrm{s}}}{2T_{\mathrm{I}}}\right)e[k-1],$$

(78)

and the step response equivalent transformation gives the difference equation of

$$u\left[k\right]=u[k-1]+K_{\mathrm{PI}}e\left[k\right]-K_{\mathrm{PI}}\left(1-\frac{T_{\mathrm{s}}}{T_{\mathrm{I}}}\right)e[k-1].$$

(79)

In each of the three specifications, the k-th value of the control signal $u\left[k\right]$ can be computed iteratively by the k-th value of the error signal $e\left[k\right]$ and the $(k-1)$-th value of the error signal and the control signal. The implementing program is, therefore, the same in all three cases, but each value must be multiplied by a different weight. For a possible Arduino code, see Appendix A, Appendix B and Appendix C.

The PI controller output signal can be approximated by

$$u\left(t\right)=K_{\mathrm{PI}}\left(e\left(t\right)+\frac{1}{T_{\mathrm{s}}}{\int }_0^{t}e\left(\tau \right)\mathrm{d}\tau \right),$$

(80)

which results in the difference equation:

$$u\left[k\right]=K_{\mathrm{PI}}\sum _{i=0}^{k-1}e\left[i\right]=u[k-1]+K_{\mathrm{PI}}e\left[k\right].$$

(81)

The control signal by the PD controller can be calculated by

$$u\left[k\right]=\frac{T_{\mathrm{D}}^{\prime }}{T_{\mathrm{s}}+T_{\mathrm{D}}^{\prime }}u[k-1]+K_{\mathrm{PD}}\frac{T_{\mathrm{s}}+T_{\mathrm{D}}}{T_{\mathrm{s}}+T_{\mathrm{D}}^{\prime }}e\left[k\right]-K_{\mathrm{PD}}\frac{T_{\mathrm{D}}}{T_{\mathrm{s}}+T_{\mathrm{D}}^{\prime }}e[k-1],$$

(82)

when applying the approximation $s=\frac{z-1}{T_{\mathrm{s}}z}$. The Tustin formula $s=\frac{2}{T_{\mathrm{s}}}\frac{z-1}{z+1}$ results in

$$u\left[k\right]=-\frac{T_{\mathrm{s}}-2T_{\mathrm{D}}^{\prime }}{T_{\mathrm{s}}+2T_{\mathrm{D}}^{\prime }}u[k-1]+K_{\mathrm{PD}}\frac{T_{\mathrm{s}}+2T_{\mathrm{D}}}{T_{\mathrm{s}}+2T_{\mathrm{D}}^{\prime }}e\left[k\right]+K_{\mathrm{PD}}\frac{T_{\mathrm{s}}-2T_{\mathrm{D}}}{T_{\mathrm{s}}+2T_{\mathrm{D}}^{\prime }}e[k-1],$$

(83)

and finally, the step response equivalent is as follows

$$u\left[k\right]=qu[k-1]+K_{\mathrm{PD}}\frac{T_{\mathrm{D}}}{T_{\mathrm{D}}^{\prime }}e\left[k\right]+K_{\mathrm{PD}}\left(1-q-\frac{T_{\mathrm{D}}}{T_{\mathrm{D}}^{\prime }}\right)e[k-1],$$

(84)

where $q={\mathrm{e}}^{-T_{\mathrm{s}}/T_{\mathrm{D}}^{\prime }}$.

Also in these specifications, the k-th value of $u\left[k\right]$ can be computed iteratively by the k-th value of the error signal $e\left[k\right]$ and the $(k-1)$-th value of the error signal and the control signal.

The PD controller output signal can be approximated by

$$u\left(t\right)=K_{\mathrm{PD}}\left(e\left(t\right)+T_{\mathrm{D}}\frac{\mathrm{d}e\left(t\right)}{\mathrm{d}t}\right)$$

(85)

which results in the difference equation:

$$u\left[k\right]=K_{\mathrm{PD}}e\left[k\right]+K_{\mathrm{PD}}{T}_{\mathrm{D}}\frac{e\left[k\right]-e[k-1]}{T_{\mathrm{s}}}.$$

(86)

The approximate formulae of the PI-PD controller used are much longer, these are not detailed here, see the result in the code (Appendix A).

Analysis of the Continuous-Time PID Controllers

The following reference signal contains three jump-like signal variations and one linearly varying segment, which allows us to study the behavior of closed-loop control systems:

$$n_{\mathrm{ref}}=\left\{\begin{array}{ccc}1000,\hfill & \mathrm{if}& t<1\mathrm{s},\hfill \\ 2000,\hfill & \mathrm{if}& 1\mathrm{s}\le t<2\mathrm{s},\hfill \\ 2500,\hfill & \mathrm{if}& 2\mathrm{s}\le t<3\mathrm{s},\hfill \\ 4750-750t,\hfill & \mathrm{if}& 3\mathrm{s}\le t<5\mathrm{s}.\hfill \end{array}\right.$$

(87)

The reference signal tracking of the P controller with phase margins ${\phi }_{\mathrm{m}}={60}^{\circ }$ and ${\phi }_{\mathrm{m}}={40}^{\circ }$ are shown in Figure 18. The measured signal (green line) is noisy, so it is easily recognizable next to the simulated noiseless curves (red curve). The reference signal is the blue line.

The measured and simulated results show excellent agreement. The gain of the controller is $K_{\mathrm{P}}=0.1746$ for the larger phase margin and $K_{\mathrm{P}}=0.3963$ for the smaller one, i.e., the loop gain $K=32.9K_{\mathrm{P}}$ is $K=5.744$ and $K=13.038$, respectively. From this, the static error can be computed according to $h=\frac{1}{1+K}$, which is $0.148$ and $0.071$, respectively. The measured and calculated static errors are equal to three decimals. It can be clearly seen that as the loop gain is increased, i.e., the phase margin is decreased, the overshoot in the transient is increased. The overshoot can be reduced by reducing the loop gain, but then the static error increases. It can be clearly seen that static error and overshoot reduction are conflicting demands. In the linearly decreasing section, the two curves are not parallel, and the error increases beyond all limits, since the type number of the control loop is $i=0$, although the error appears to decrease (over the time interval considered) as the loop gain increases. It can be concluded that the problem cannot be solved satisfactorily with a purely proportional controller.

The reference tracking of the PI controller at ${\phi }_{\mathrm{m}}={60}^{\circ }$ is shown in Figure 19, where the reference, measured and simulated signals are the blue, green and red lines, respectively. Here, the measured and simulated results show excellent agreement. In transient, there is a small overshoot which can be reduced by decreasing the gain $K_{\mathrm{PI}}$, i.e., by increasing the phase margin, in this problem the overshoot at ${\phi }_{\mathrm{m}}={72}^{\circ }$ phase margin practically disappears (then $K_{\mathrm{PI}}=0.0431$). In the stationary state, the error is zero, since the integrating term increases the type number to $i=1$, and the speed follows the reference with a constant error in the linearly decreasing phase. It can be said that the PI controller is a good choice to solve this problem.

The reference tracking of the PD controller at ${\phi }_{\mathrm{m}}={60}^{\circ }$ is shown in Figure 20, with good agreement between the measured and simulated results. Again, reference, measured and simulated signals are the blue, green and red lines, respectively. The steady-state behavior of the control loop is inherently the same as that of the P-controller. The system cannot be accelerated significantly due to the saturation of the feedback loop, the feedback signal of the designed controller is sometimes many times higher than the upper limit of the PWM excitation signal 250, which means that acceleration does not occur in the real system. This is well observed for the time step jumps $t=0\mathrm{s}$ and $t=1\mathrm{s}$. Here, the simulated speed has an overshoot, while the measured speed is without overshoot and the run-up is below the theoretical result. The effect of the differentiating branch in the controller causes the control signal to respond to the change in the error signal over a wide range with an abrupt, unpreferred change jerking the amplifier input. The PD controller is not recommended for this application for the above reasons.

The result of the PI-PD controller is not discussed separately because it is very similar to the result of the PI controller.

To set the sampling period $T_{\mathrm{s}}$, the following rule of thumb can be used:

$$T_{\mathrm{s}}=\frac{2\pi }{180{\omega }_{\mathrm{c}}}\Delta {\phi }_{\mathrm{m}},$$

(88)

where the maximum value of the phase decay $\Delta {\phi }_{\mathrm{m}}$ is set to $5^{\circ }$. The sampling period time is $6.9\mathrm{ms}$, $11.5\mathrm{ms}$, $4.1\mathrm{ms}$, $4.1\mathrm{ms}$ and $1.9\mathrm{ms}$ for P, PI, PD, PI-PD ($n=2$), and PI-PD ($n=5$) controllers according to the cut-off frequency data of Table 2, respectively. By simply running the microcontroller’s $\mathtt{loop}\left(\right)$ function, the resulting sampling period time for this implementation is on average around 30 ms, which corresponds to a sampling frequency of roughly 30 Hz, and is larger than the values given by the rule of thumb. However, it can be stated that the controllers perform satisfactorily, very close to the theoretical results. The largest difference in the calculated and measured transient was observed for the PD controller.

Design of the Discrete-Time PID Controllers

When designing a discrete-time PID controller, the first step is to transform the continuous-time transfer function of the plant into a step response equivalent form, where the sampling period $T_{\mathrm{s}}$ is predefined. The discrete-time transfer function of a DC motor is generally as follows:

$$W_{\mathrm{P}}\left(z\right)=M\frac{z-z_1}{(z-p_1)(z-p_2)},$$

(89)

where M is a constant, and the discrete-time model has one zero $z_1$ and two poles $p_1$ and $p_2$. These of course depend on the sampling period. For the transformation, the Matlab R2015a function $\mathtt{c}\mathtt{2}\mathtt{d}$ is used. As an example, the transfer function at 50 Hz and 100 Hz is as follows:

$$\begin{array}{cc}\hfill W_{\mathrm{P}}\left(z\right)=& 0.8805\frac{z+0.8055}{(z-0.8817)(z-0.5922)},\hfill \\ \hfill W_{\mathrm{P}}\left(z\right)=& 0.2442\frac{z+0.8972}{(z-0.9390)(z-0.7695)}.\hfill \end{array}$$

(90)

The poles of a discrete-time system can be mapped from the poles of the continuous-time system according to the relation $p_{\mathrm{discrete}}={\mathrm{e}}^{p_{\mathrm{continuous}}{T}_{\mathrm{s}}}$.

The design is conducted on a phase margin, in a similar way to the design of continuous-time controllers, by solving the following system of nonlinear equations:

$$\left|W_{\mathrm{O}}\left({\mathrm{e}}^{\mathrm{j}{\omega }_{\mathrm{c}}{T}_{\mathrm{s}}}\right)\right|=1,-{180}^{\circ }+{\phi }_{\mathrm{m}}=\mathrm{arc}{W}_{\mathrm{O}}\left({\mathrm{e}}^{\mathrm{j}{\omega }_{\mathrm{c}}{T}_{\mathrm{s}}}\right).$$

(91)

The system of equations is solved by calling the Matlab R2015a function $\mathtt{fsolve}$. The step response of the designed closed-loop control systems are compared in Figure 21, and some characteristic data are summarized in

Table 3. Some data of closed-loop systems for the phase margin ${\phi }_{\mathrm{m}}={60}^{\circ }$, the sampling frequency is 100 Hz.

Controller	${\mathit{\omega }}_{\mathbf{c}}[\mathbf{rad}/\mathbf{s}]$	K	${\mathit{v}}_{max}$	${\mathit{v}}_{\mathbf{st}}$	$\mathbf{\Delta }\mathit{v}[\%]$	${\mathit{u}}_{\mathbf{t}}=\frac{\mathit{u}\left(0\right)}{{\mathit{u}}_{\mathbf{st}}}$	${\mathit{u}}_{\mathbf{st}}$
P	21.84	0.143	0.98	0.83	18.9	5.78	0.025
PI	12.96	0.072	1.08	1.00	8.89	2.49	0.031
PD	63.66	1.385	1.05	0.91	14.6	49.9	0.027
PI-PD	53.85	1.191	1.07	1.00	7.28	39.2	0.031

. The nature of the step responses is of course the same as the results obtained for the continuous-time design (see Figure 17).

By reducing the sampling period time, the results of the discrete-time model converge to the results of the continuous-time model, while the value of the control signal increases at the beginning of the transients.

The result of the controller design is a difference equation, which can be derived, for example, for the PI controller as follows. Start from the transfer function of the controller:

$$W_{\mathrm{C}}\left(z\right)=\frac{U\left(z\right)}{E\left(z\right)}=K_{\mathrm{PI}}\frac{z-z_{\mathrm{I}}}{z-1}=K_{\mathrm{PI}}\frac{1-z_{\mathrm{I}}{z}^{-1}}{1-z^{-1}},$$

(92)

from which, after cross multiplication, the following equation can be written:

$$K_{\mathrm{PI}}E\left(z\right)-K_{\mathrm{PI}}{z}_{\mathrm{I}}E\left(z\right)z^{-1}=U\left(z\right)-U\left(z\right)z^{-1},$$

(93)

to which the following difference equation is given by the shift operator of the z-transform:

$$u\left[k\right]=u[k-1]+K_{\mathrm{PI}}e\left[k\right]-K_{\mathrm{PI}}{z}_{\mathrm{I}}e[k-1].$$

(94)

This equation gives the algorithm of the controller.

The PD controller can be implemented in a similar way as above:

$$u\left[k\right]=K_{\mathrm{PD}}e\left[k\right]-K_{\mathrm{PD}}{z}_{\mathrm{D}}e[k-1],$$

(95)

and the algorithm of the PID controller is as follows:

$$u\left[k\right]=u[k-1]+K_{\mathrm{PID}}e\left[k\right]-K_{\mathrm{PID}}(z_{\mathrm{I}}+z_{\mathrm{D}})e[k-1]+K_{\mathrm{PID}}{z}_{\mathrm{I}}{z}_{\mathrm{D}}e[k-2].$$

(96)

Thus, the k-th value of the control signal is the weighted sum of the $(k-1)$-th value of the control signal, the k-th, $(k-1)$-th and $(k-2)$-th values of the error signal. The controller program is very easy to implement using a general form of

$$W_{\mathrm{C}}\left(z\right)=\frac{b_0{z}^2+b_{1}z+b_2}{z^2+a_{1}z+a_2}=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1{z}^{-1}+a_2{z}^{-2}},$$

(97)

which results in the difference equation:

$$u\left[k\right]=-a_{1}u[k-1]-a_2[k-2]+b_{0}e\left[k\right]+b_{1}e[k-1]+b_{2}e[k-2].$$

(98)

The included coefficients are summarized in

Table 4. Parameters of the difference Equation (98), ${\phi }_{\mathrm{m}}={60}^{\circ }$ and $f_{\mathrm{s}}=50\mathrm{Hz}$.

Controller	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{b}}_0$	${\mathit{b}}_1$	${\mathit{b}}_2$
P	0	0	0.1241	0	0
PI	1.00	0	0.0642	−0.0566	0
PD	0	0	0.4568	−0.2705	0
PI-PD	1.00	0	0.3577	−0.5272	0.1868

and in

Table 5. Parameters of the difference Equation (98), ${\phi }_{\mathrm{m}}={60}^{\circ }$ and $f_{\mathrm{s}}=100\mathrm{Hz}$.

Controller	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{b}}_0$	${\mathit{b}}_1$	${\mathit{b}}_2$
P	0	0	0.1430	0	0
PI	1.00	0	0.0720	−0.0676	0
PD	0	0	1.3851	−1.0659	0
PI-PD	1.00	0	1.1911	−2.0351	0.8607

at 50 Hz and 100 Hz sampling frequencies, respectively.

The results for discrete-time PID controllers are not discussed separately but are very similar to those for continuous-time PID controllers. The  implemented Aruino code of $f_{\mathrm{s}}=50\mathrm{Hz}$ is presented in Appendix B for further use.

3.3.2. State Feedback Controllers

Observer Design to Measure the Current Signal

The discrete-time state-space model of the continuous-time system is obtained by the Matlab R2015a function $\mathtt{c}\mathtt{2}\mathtt{d}$. The sampling period is defined to $T_{\mathrm{s}}=1/50\mathrm{s}$ ($f_{\mathrm{s}}=50\mathrm{Hz}$), and 

$$\mathbf{x}[k+1]=\left[\begin{array}{cc}0.452& -0.001956\\ 13.9& 0.9552\end{array}\right]\mathbf{x}\left[k\right]+\left[\begin{array}{c}0.08457\\ 0.9606\end{array}\right]u\left[k\right],y\left[k\right]=\left[\begin{array}{cc}0& 1\end{array}\right]\mathbf{x}\left[k\right].$$

(99)

The eigenvalues of the system matrix, i.e., the poles of the system, are $p_1=0.8937$ and $p_2=0.5135$.

The vector $\mathbf{G}$ is determined by Ackermann’s formula, with the observer’s poles at $0.05$. The Matlab R2015a function $\mathtt{acker}$ gives the result of

$$\mathbf{G}=\left[\begin{array}{c}0.0257\\ 0.9946\end{array}\right],\mathbf{H}=\left[\begin{array}{c}0.0846\\ 0.9606\end{array}\right],\mathbf{F}=\left[\begin{array}{cc}0.4520& -0.0277\\ 13.8997& -0.0393\end{array}\right].$$

(100)

The approximation of the armature current and speed for a stepped excitation is shown in Figure 22. The measured current signal is quite noisy in this range, but the currents calculated with the continuous-time observer and the discrete-time observer are close to each other. The measured and observed speeds are in perfect agreement.

Kalman Filter to Estimate the Current Signal

To determine the covariance matrix of the Kalman filter ${\mathbf{R}}_{\mathrm{v}}$ and the covariance $R_{\mathrm{w}}$, the measurement result shown in Figure 23 is used when the motor is excited with a constant PWM signal of 150.

The mean values calculated from the $n=201$ sample are ${\mu }_{i_{\mathrm{a}}}=0.155$, ${\mu }_{\omega }=222.2334$, and the covariances are as follows:

$${\mathbf{R}}_{\mathrm{v}}=\left[\begin{array}{cc}0.0266& -0.0209\\ -0.0209& 2.5720\end{array}\right],R_{\mathrm{w}}=2.572.$$

(101)

The gain of the Kalman filter $\mathbf{G}$ can be determined by the Matlab R2015a function $\mathtt{dare}$, which numerically solves the Riccati Equation (44). Then, $\mathbf{F}$ and $\mathbf{H}$ can be calculated:

$$\mathbf{G}=\left[\begin{array}{c}0.0043\\ 0.9583\end{array}\right],\mathbf{H}=\left[\begin{array}{c}0.0846\\ 0.9606\end{array}\right],\mathbf{F}=\left[\begin{array}{cc}0.4520& -0.0062\\ 13.8997& -0.0031\end{array}\right].$$

(102)

The estimated speed signal and the estimated current signal obtained by the Kalman filter have been used in the control loop.

Design of State Feedback with Pole Placement and Ackermann’s Formula

The eigenvalues of the extended system matrix, i.e., the poles of the extended system, are $p_1=0.8937$, $p_2=0.5135$, and $p_3=1$. The last one is due to the integrating effect, the other two poles are from the discrete-time model of the motor.

Two cases are described in more detail.

In the first case, the poles $p_1$ and $p_2$ are not moved, the pole $p_3=1$ is changed to $p_3=1/2$, i.e., the pole connected to the integrator is accelerated. This gives the following gains using Ackermann’s formula (see Figure 24, where reference, measured and simulated signals are the blue, green and red lines):

$${\mathbf{p}}^{\mathrm{T}}=\left[\begin{array}{ccc}{p}_1& p_2& p_3/2\end{array}\right]\Rightarrow {\mathbf{k}}^{\mathrm{T}}=\left[\begin{array}{ccc}4.0835& 0.1610& 0.7598\end{array}\right].$$

(103)

As a result of the experimentation with pole placement, further modification of $p_3$ (e.g., to $1/4$ or $1/5$) has practically no effect on the behavior of the closed loop. Changing the pole $p_2$ also has no effect on the controlled characteristic or the control signal. Moving the pole $p_1$ is more significant. At the following values, as in the second case, the control signal does not exceed the upper limit, but a faster settling is obtained compared to the previous one (see Figure 25, where reference, measured and simulated signals are the blue, green and red lines):

$${\mathbf{p}}^{\mathrm{T}}=\left[\begin{array}{ccc}{p}_1/1,5& p_2& p_3/2\end{array}\right]\Rightarrow {\mathbf{k}}^{\mathrm{T}}=\left[\begin{array}{ccc}5.9413& 0.3076& 2.8886\end{array}\right].$$

(104)

Further acceleration, however, will result in overshooting or saturation of the actuating block.

LQ-Based State Feedback Design

After some trials, the following choice leads to a very similar result to the second case obtained by Ackermann’s formula:

$$\mathbf{Q}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 100\end{array}\right],r=100\Rightarrow {\mathbf{k}}^{\mathrm{T}}=\left[\begin{array}{ccc}2.6691& 0.1270& 0.8351\end{array}\right].$$

(105)

Controller gains are obtained by the solution of (52) and (53).

A possible discrete-time state feedback controller implementation can be studied in Appendix C for further use.

4. Discussion

The DC motor analysis with real laboratory tests showed the well-known theoretical result that the DC motor can be modeled as a linear plant in the control systems. From the steady-state and switch-on measurements, the linear state space model as well as the transfer function model have been identified, and it is supposed that the model is second order. Measurement results and simulated curves show a very good agreement.

Continuous-time PID controllers have been designed via the well-known phase margin method, and the discretized form of the controllers has been implemented in the Arduino UNO microcontroller. The discrete-time model-based PID controllers are also designed and realized. The two approaches result in the same control signal as well as the same output signal. The continuous-time model is a little bit easier to understand, and the results of the discrete-time model proved that the continuous-time controller can be used efficiently; however, the discrete-time model is closer to the theory of discrete-time controller design. Moreover, discrete-time modeling can be implemented in microcontrollers. More complicated problems should be worked out in the frame of the discrete-time theory because the sampling time can be fixed at the very first step of the design.

The noisy current signal is efficiently measured by the state observer as well as by the Kalman filter. It is highlighted that the Kalman filter did not result in a much better noise reduction, meaning that the Kalman filter parameters should be adjusted as a future task as well as the noise model parameters. Practically, in this application the deterministic state observer and the Kalman filter gave the same measurement results.

State feedback controller design has been realized and the pole placement has been studied comprehensively. LQ optimal controller was also studied. The design of the parameters, i.e., the matrix $\mathbf{Q}$ and the parameter R should be improved, to reach a predefined closed loop system output and control signal [55].

There are some other linear controllers that can be realized, e.g., the two degrees of freedom controller, dead-beat controller, optimal controller with conditions on the control signal, model predictive controller and so forth. Linear control suffers from obvious limitations in terms of lack of disturbance handling and modeling uncertainties. The effect of uncertain model parameters and disturbances can be decreased by the integral term of the control loops mentioned in the paper; however, uncertainty and disturbance can be handled efficiently by robust control methods [56,57]. This paper has not worked on these questions, these can be the focus of fruitful future tasks. The paper mentioned only a full-order observer to estimate the state vector elements; however, reduced-order state estimation can be applied. This can be a future task as well.

The test bench is easy to realize, the motor, the sensors and the microcontroller are very easy to build and assemble, and it can be a useful tool for engineering students in undergraduate systems and control lectures and laboratories. This is why the Arduino codes are inserted for further studies.

There are, of course, many other microcontroller-based platforms for solving the problem presented in this review. The Arduino UNO is a cheap and simple solution that is also very useful in undergraduate control engineering education. For more complex problems that require more computing power or special tasks, other types of hardware should be used.