*Proceeding Paper* **Automating the Study of a Linearized Model of Diabetes Mellitus and Tuning a PID Controller for This Model †**

**Denis Andrikov and Sinan Kurbanov \***

The Patrice Lumumba Peoples' Friendship University of Russia, 6 Miklukho-Maklaya St., Moscow 117198, Russia; andrikov-da@rudn.ru

	- 14–16 December 2022.

**Abstract:** This paper proposes a simulation linear model created in the MATLAB environment, which provides a process for regulating blood sugar levels. The controller is built for the need for any type of diabetes to control and normalize the blood sugar content of the patient in order to eliminate the differences in the quality of life of a diabetic patient and a healthy person. The linearization of the nonlinear model was performed, and the adequacy of the linearized model was verified and confirmed using the MATLAB simulation. The choice of the PID controller and the CHR method for its adjustment was justified and MATLAB tools were used to show the implementation of these methods. The model of the patient with the controller has been built; the algorithm for the automatic adjustment of the PID controller parameters has been developed and realized. The directions for continuation of the work on this problem regarding regulation in the system under study are proposed.

**Keywords:** diabetes mellitus; insulin pump; automatic control system; insulin–glucose interaction models; PID controller; linearization; parameter setting

### **1. Introduction**

Diabetes mellitus (DM) is a group of metabolic diseases characterized by chronic hyperglycemia, which is the result of impaired insulin secretion, the action of insulin, or both. Chronic hyperglycemia leads to damage and dysfunction of various organs, including the eyes, kidneys, nerves, heart, and blood vessels. In any type of diabetes, the control and normalization of blood sugar becomes one of the main tasks of the patient and his attending physician. When the sugar level is closer to the normal limits, the risk of complications is lower, and there are fewer symptoms of diabetes as well as fewer differences between the quality of life of a patient with diabetes and a healthy person. The current state of medicine and science in general allows the treatment of diabetes mellitus on an outpatient basis, mainly through the use of an individual control system, that is, an insulin pump. The algorithms of such a system can be implemented using a controller.

According to research published in the tenth edition of the Diabetes Atlas of the International Diabetes Federation (IDF), 537 million people aged from 20 to 79 years have diabetes, and this number is steadily increasing every year [1], including in the Russian Federation [2]. This is why the topic of developing new and improving existing automatic control systems for regulating blood sugar levels in patients with diabetes mellitus is becoming increasingly relevant. The aim of the work is the practical implementation of a simulation model in MATLAB for regulating blood sugar levels using an actuator, namely an insulin pump, in order to be able to use this technique when such compact medical devices are put into operation. To achieve this goal, the following tasks were set and solved:

1. Study and analysis of existing mathematical models of diabetes mellitus in terms of applying the methods of the theory of linear automatic control systems to them [3–15].

**Citation:** Andrikov, D.; Kurbanov, S. Automating the Study of a Linearized Model of Diabetes Mellitus and Tuning a PID Controller for This Model. *Eng. Proc.* **2023**, *33*, 42. https://doi.org/10.3390/ engproc2023033042

Academic Editors: Askhat Diveev, Ivan Zelinka, Arutun Avetisyan and Alexander Ilin

Published: 27 June 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).


#### **2. Materials and Methods**

#### *Rationale for Choosing a Mathematical Model*

After the possibility of an external control for blood glucose levels in patients with type I diabetes mellitus by monitoring blood glucose levels and injecting insulin was experimentally proven, it became necessary to develop a safe and effective algorithm for controlling the devices that solve the problem of diabetes for patients. When developing a mathematical model of the interaction between insulin and glucose in the human body, pharmacokinetic and physiological models were considered. The main difference between pharmacokinetic models and physiological ones is that they contain fewer parameters, so they are easier to use for interpreting experimental data. Physiological models are characterized, as a rule, by a large number of parameters. Therefore, Sorensen's model [8] has 21 states and 22 metabolic functions that describe the dynamics of glucose, insulin, and glucagon. This model is widely used in glucose monitoring and has been criticized for inaccurately representing the observed changes in glucose [9]. Another physiological model, from Kobelli [10], is nonlinear and consists of glucose, insulin, and glucogenic subsystems; it has nine states, 23 metabolic functions, and 46 parameters. Therefore, only some pharmacokinetic models were considered as prototypes for modeling: the Bergman model [11], the Hovorka model [12], and the simplest basic differential mathematical model. The simplest basic differential mathematical model is of the following form:

$$\begin{cases} \frac{dx}{dt} = -a\_1 \mathbf{x}y + a\_2(\mathbf{x}\_0 - \mathbf{x}) + a\_3 P(t) \\ \frac{dy}{dt} = b\_1(\mathbf{x} - \mathbf{x}\_0)H(\mathbf{x} - \mathbf{x}\_0) - b\_2 y + b\_3 u(t) \end{cases} \tag{1}$$

where *x*(*t*) is the blood sugar level, *y*(*t*) is the level of insulin in the blood, and *x*<sup>0</sup> is the fasting blood sugar level. The constants *a*1, *a*2, and *a*3, as well as *b*1, *b*2, and *b*3, are positive and are the sensitivities of the sugar and insulin gradients, respectively. The function *H*(*x*) is a single step, *P*(*t*) describes changes in sugar levels from food intake, and *u*(*t*) describes changes in insulin levels. The undoubted advantages of this model are its simplicity and compactness, consisting in a minimum number of equations and parameters. The adequacy of this model is proved in practice [14]. The model has two input variables, *P*(*t*) and *u*(*t*), and two state variables, *x*(*t*) and *y*(*t*). The state variables are the main variables of the model. They are quantities that uniquely determine the current state of the control object. These variables can be changed (i.e., controlled) in clinical practice. As follows from the equations in (1), the mathematical model is nonlinear. The equations in (1) include the product of the state variables; in both equations, the coefficients on the state variable x depend on this variable, and this dependence is also nonlinear because one of the comultipliers is a step function. The deviations of the state variables from the initial value ("starvation") are commensurate with the value itself [14,15] and, therefore, cannot be considered small. In this connection, the decomposition of nonlinear functions into a Taylor series for the purpose of linearization of the object will be ineffective because, to achieve acceptable accuracy, the account of a large number of row members will be required. The attempt to linearize differential equations at steady-state constant values for "hunger" *x*0, *y*<sup>0</sup> via a Taylor series expansion is impossible because of the infinite derivatives of the step functions *H*(*x* − *x*0) at point *x*0. As stated in [16], for this kind of linearization, the derivatives must have a single and finite value; otherwise, the equation is nonlinearizable. The application of I.A. Vyshnegradsky's hypothesis for small deviations of all variables from steady-state values seems to be an acceptable method of linearization [16]. The steadystate value is logically considered to be the function of the state variables for a healthy person. The model was linearized by entering a new variable—the deviation of the sugar level from zero, corresponding to hunger:

$$z = x - x\_0$$

and the following matrices of the linear system in the state space were obtained

$$\mathbf{A} = \begin{pmatrix} -a\_2 H(-z) & -a\_1 \mathbf{x}\_0 \\ b\_1 H(-z) & b\_2 \end{pmatrix}; \mathbf{B} = \begin{pmatrix} a\_3 & 0 \\ 0 & b\_3 \end{pmatrix}; \mathbf{C} = \begin{pmatrix} 1 & 0 \end{pmatrix}; \mathbf{D} = \begin{pmatrix} \mathbf{x}\_0 \end{pmatrix}.$$

The transfer function with a second-order characteristic polynomial was derived from the representation in the state space

$$Z(s) = -\frac{a\_1 \ge\_0 b\_3}{s^2 + b\_2s + a\_1 \ge\_0 b\_1} \mathcal{U}(s) + a\_3 \frac{s + b\_2}{s^2 + b\_2s + a\_1 \ge\_0 b\_1} P(s)$$

The simulation showed a practical coincidence of the transition functions of the linearized and the original system within two to three hours from the moment of eating (Figure 1), so, the linear model was accepted for further research provided that the maximum of the transition function of the system with the controller is reached at the specified time.

**Figure 1.** Comparison of transients in the initial (*z*,*y*) and linearized (*z*1,*y*1) models. The control object is a patient with diabetes mellitus (DM) [17].

When choosing the controller configuration method, the main preference was given to engineering methods that allow for calculation of the controller parameters based on the characteristics of the transient process graph. To date, the most widespread use is made of PID controllers. Practical experience shows that the use of PID controllers in complex systems (devices) provides a sufficiently small control error and the transition process required by the operating conditions of most control objects. Therefore, the PID regulation will be the priority in the control task under study. Currently, there are quite a lot of different methods for tuning the parameters of PID controllers; still, the most common is the Ziegler–Nichols tuning method [18], proposed by its authors in 1942. This method is simple in application but gives not very good results. As a rule, after setting the regulator parameters, one has to make manual adjustments in order to improve the quality of regulation. Nevertheless, this method is still often used in practice, although many more accurate methods have become available. Ziegler and Nichols proposed two methods for tuning PID regulators [19]. One of them is based on the response parameters of the object to a single jump, and the other is based on the frequency characteristics of the control object. To calculate the PID controller parameters according to the first Ziegler–Nichols method, only two parameters are used: the negative coordinate "a" of the point of intersection of the tangent to the transient curve with the maximum slope, and the coordinate "L" of the point of intersection of the same tangent and the time axis. The formulas for calculating PID controller coefficients are summarized in Table 1.


**Table 1.** Formulas for calculating regulator coefficients by the Ziegler–Nichols method.

The Ziegler–Nichols method does not take into account in any way the requirements for the stability margin of the system, which is its significant disadvantage. In addition, the Ziegler–Nichols method gives parameters that are far from optimal. The second Ziegler– Nichols method (the frequency method), as input data for calculations, uses the frequency *ω*180, at which the phase shift in the open circuit reaches 180°, and the modulus of the loop gain at this frequency is *K*180. Knowing the parameter *ω*180, first, find the period of natural oscillations of the system, and then, from Table 1, determine the parameters of the regulator. The accuracy of the controller tuning and the disadvantages of both Ziegler–Nichols methods are the same. The CHR method, which uses a similar method for calculating the parameters of the PID regulator, allows us to bypass the drawbacks of the Ziegler–Nichols method. The authors of this method, Chien, Hrones, and Reswick (CHR) [20], in contrast to Ziegler and Nichols, used the criterion for the maximum slew rate in the absence of overshoot or with no more than 20% overshoot. This criterion gives a larger stability margin than in the Ziegler–Nichols method. The CHR method gives two different systems of regulator parameters. One is obtained by observing the response to the setpoint change (Table 2), and the other is obtained by observing the response to external disturbances (Table 3). Which parameter system to choose depends on what is more important for a particular controller: the quality of control when changing the setpoint or the attenuation of external disturbances. If both are important, it is necessary to use regulators with two degrees of freedom [21].

**Table 2.** Formulas for calculating regulator coefficients by CHR method, by response to set point change.



**Table 3.** Formulas for calculating regulator coefficients by CHR method, by response to external perturbations.

The CHR method uses an approximation of the object by a first-order delayed model. CHR uses the same initial parameters a and L as the Ziegler–Nichols method. Thus, each system requires a different approach to adjust the controller parameters. A comparative analysis of the methods was carried out, and the CHR method was adopted as providing a better quality of the transient process and a greater margin of stability of the system with a controller compared to the Ziegler–Nichols method. A method was developed for determining the parameters of the PID controller based on a single input variable—the angle of inclination of the tangent to the transient graph. The method was based on the approximation of the transition process graph by a polynomial, followed by analytical differentiation of this polynomial, in order to accurately determine the angle of the tangent. The following algorithm of PID controller tuning was developed:


Step 2 of the algorithm—the task of tangent construction, or the task of numerical differentiation of a tabularly defined function—requires detailed consideration. The simplest solution to this problem is based on determining the derivative of the function:

$$y' = \lim\_{\Delta x \to 0} \frac{\Delta y}{\Delta x}$$

and uses various finite difference formulas [22]; however, the finite difference method is associated with known difficulties in estimating the error of the result. As stated in [22], the main components of the error of numerical differentiation are the approximation error (also called the truncation error) and rounding errors in computer calculations. The approximation error is determined by the magnitude of the residual term—the difference between the approximated and actual values of the derivative. It is noted in [22] that the analysis of the residual term is nontrivial and that the approximation error tends to decrease as the step Δ*x* decreases. In contrast to the approximation error, the rounding error increases as the step Δ*x* decreases. Therefore, the total error of numerical differentiation can decrease as the step decreases only up to some limiting value, after which a further decrease in the step will, at best, not improve the accuracy of the results. In addition, finite difference formulas usually use a constant step, and the resulting simulation function may be a variable step. Optimal accuracy can be achieved by regularizing the numerical differentiation procedure. The simplest way to regularize is to choose a step Δ*x* such that the inequality is as follows:

$$|f(\mathbf{x} + h) - f(\mathbf{x})| > \varepsilon,\tag{2}$$

where *ε* > 0 is some small number. When calculating the derivative, this eliminates the subtraction of closely related numbers, which usually leads to an increase in errors. This is all the more dangerous when subsequently dividing the increment of a function by a small number Δ*x*. This method cannot be applied to our problem for the reason mentioned above for the variable step Δ*x* in the original table function.

Another method for regularization is to smooth the tabulated values of the function by selecting some smooth approximating function, for example, a polynomial. Further differentiation of the polynomial using the derivatives table does not present computational difficulties and does not introduce additional errors. This method was used to solve the problem.

The polynomial calculated by the built-in polyfit function will be used as the interpolation polynomial. This choice among many variants [22]—Lagrange, Hermite, and other polynomials—is dictated by the simplest implementation in the program. The transient graph is a fairly smooth curve. By plotting the polynomial trend on this graph, the degree of the polynomial that showed the best coincidence with the graph of the original function was determined. It was found that the polynomial should be of degree 10. This does not exclude a lower degree polynomial; the change in degrees was implemented in the program. The following algorithm was developed to form a set of graph points to build the approximating polynomial:


After forming a vector of approximating polynomial coefficients by this set, this vector was transformed into a vector of derivative coefficients according to the following scheme:

$$
\begin{pmatrix}
polynomial \\ a\_n \\ a\_{n-1} \\ \dots \\ a\_1 \\ a\_0 \end{pmatrix} \to \begin{pmatrix}
 derivative \\ na\_n \\ (n-1)a\_{n-1} \\ \dots \\ 1a\_1 \\ 0 \end{pmatrix} \to \begin{pmatrix}
 derivative \\ na\_n \\ (n-1)a\_{n-1} \\ \dots \\ 1a\_1 \end{pmatrix}
$$

and the value of this derivative *y* (0) at the point of food intake point was calculated. Since the food intake point has the coordinates (0;0), the equations of the tangent are as follows:

$$y(x) = y'(0) - x$$

and

$$a = y'(0) - L$$

The parameter *L* has a clear biophysical meaning and is the average time elapsed from the beginning of a meal to the feeling of fullness, which is the signal of an increase in blood sugar level because of the digestion of nutrients. The study of this parameter and the scientific substantiation of its value is a topic for a separate study that is beyond the scope of this work, so we took an a priori value of 15 min or

$$L = 0.25 \text{h.}$$

Each control system requires its own approach for setting the controller parameters. In this case, the setup process is hampered by some factors, namely:


Considering these factors, the MATLAB/Simulink software package was selected to study the linearized insulin injection control system and test the PID controller tuning algorithm. The development of the model in Simulink is characterized by relative simplicity and clarity, and the management of this model from the Matlab script (m-file) allows for automation of the process of setting model parameters and processing simulation results. In the process of research, it was necessary to take into account a significant variation in the transfer function coefficients of the model. Therefore, the model of diabetes mellitus itself was selected in the mathematical modeling program. Its visual implementation was carried out by means of Simulink. The model is shown in Figure 2. In accordance with the type of transfer function, the model has two channels—one for nutrition and another for insulin injection. The transfer function of each channel is recorded with parametrization through the components of the coefficient vectors *a*(*i*) and *b*(*i*). The PID controller as a block from the Simulink library is included in the feedback loop of the channel for insulin injection. Physically, the PID controller is an insulin pump. The nutrition function according to [4] has the form

$$P(t) = \begin{cases} 0, & t < t\_0 \\ Q e^{-K(t - t\_0)}, & t \ge t\_0 \end{cases}$$

where *Q* is the amount of food, *K* is a parameter that characterizes the type of sugar that comes with food, and *t*<sup>0</sup> is the meal time. This function is implemented in the MATLAB function block. The signals from the channel outputs are summed according to the expression for the transfer function *Z*(*s*). The transfer to MATLAB of the following simulation results is provided—the model time is represented as *t*, and the insulin level is represented as *y*.

**Figure 2.** Development of a control object model in MATLAB/Simulink [17].

The interaction algorithm between the model and the MATLAB script that calls it includes the following steps:


Then, steps 1–3 are repeated according to the purpose of the study; i.e., these items are executed the required number of times or until another termination condition is met [17].

#### **3. Results & Discussions**

The criterion for completing the adjustment was the specified excess of the maximum insulin level over the level of a healthy person, while the level of a healthy person was calculated by a direct numerical solution of the original nonlinear differential in Equation (1). A number of computational experiments were carried out on the model, as a result of which the limits of changing the parameters of the PID controller were determined:


A transient graph with the controller settings *K* = 18, *Ti* = 0.0333 h, and *Td* = 0.0058 h, for a patient with an insulin pump (IP), is shown in Figure 3. For comparison, the same graph shows the transition process in the "healthy" model.

As a result of the work completed, the prospects of applying the considered approach to the development of controllers for insulin pumps can be considered confirmed. This is evidenced by the results of modeling and comparing the parameters of the transition process (see Figure 3). The advantage of using a linear model is the simplification of the controller tuning procedure, and linearization allows you to automate this procedure. Thus, it is possible to achieve the dynamics of the blood sugar level of the "patient" to those characteristic of the "healthy".

**Figure 3.** Demonstration of transient modeling results [17].

#### **4. Conclusions**

To sum up, we can say that the goal of regulation set in the work has been achieved—a simulation model has been created in MATLAB as a result of the following stages:

1. Based on the study and analysis of existing mathematical models of diabetes mellitus, in terms of the application (to them) of methods of the theory of linear automatic

control systems, a non-linear model was chosen that is consistent with the known physiological facts while having a minimal set of equations and parameters.


These results suggest the following directions for further research:


**Author Contributions:** Conceptualization, D.A.; Formal analysis, S.K.; Investigation, S.K.; Methodology, D.A.; Project administration, D.A.; Software, S.K.; Validation, D.A.; Visualization, S.K.; Writing original draft, S.K.; Writing—review and editing, D.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of The Patrice Lumumba Peoples' Friendship University of Russia (protocol code 200, 2006).

**Informed Consent Statement:** Written informed consent has been obtained from the patient(s) to publish this paper.

**Data Availability Statement:** The study did not report any data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
