A robust fixed point transformation-based approach for type 1 diabetes control

Modeling and control of diabetes mellitus (DM) are difficult due to the highly nonlinear attitude, time-delay effects, the impulse kind input signals and the lack of continuously available blood glucose (BG) level to be regulated. Regarding the mentioned problems, identification of DM model is crucial. Furthermore, due to the lack of information about the internal states (which cannot be measured in everyday life) and because the BG level is not available in every moment over time, adaptive robust control design method regardless exact model dependency would successfully handle these unfavorable effects without simplifications. The recently developed nonlinear robust fixed point transformation (RFPT)-based controller design method requires only a roughly approximate model in order to realize the controller structure. Moreover, parallel simulated approximate models—in order to provide additional internal information—can be used with the method. In this paper, the usability of the novel RFPT-based technique is demonstrated on the physiological problem of diabetes.


Introduction
The concept of modern control technology, under the name "Theory of Governors," originates from a paper by Maxwell [25], in which, without entering into any details of the particular mechanisms that were known at that time, he directed the attention of engineers and mathematicians to a more general dynamical theory. Following his fundamental achievements in the description of the electrical and magnetic phenomena [26], by the use of the elementary circuit components as resistors, capacitors and inductances as linear timeinvariant (LTI) elements, in the field of electrical engineering rapid development was produced that intensively utilized the mathematical achievements of the nineteenth century as the analysis of complex numbers, Laplace, Z, Fourier, Mellin and other transforms [6]. In this great flourishment of linear control technology, the use of, and thinking on the basis of the frequency picture became prevailing. This general attitude lasted till the beginning of 1960s when, according to [19], Rudolf Kalman "…challenged the accepted approach to control theory of that period, limited to the use of Laplace transforms and the frequency domain, by showing that the basic control problems could be studied effectively through the notion of the state of the system that evolves in time according to ordinary differential equations in which control appears as parameters. …Liberated from the confines of the frequency domain and further inspired by the development of computers, automatic control theory became the subject matter of a new science called systems theory." This liberation from the LTI systems tailored frequency domain-based problem tackling that opened the way for studying nonlinear dynamical systems in "system theory," happened relatively lately in comparison with other revelations of fundamental significance in life sciences and chemistry, regarding biophysical and/or biochemical systems.
Lapicque [22] elaborated a strongly nonlinear model for describing neuron spiking. Hodgkin and Huxley [15] quantitatively modeled the membrane current in nerve excitation. In the early 1960s, further research results were published for the pulse transmission and membrane models [10,27]. To ease the study of these nonlinear dynamic phenomena, L. Chua and T. Matsumoto constructed and studied a special nonlinear, but relatively simple electrical circuit [24]. In the last decade of the twentieth century, the significance of chaotic phenomena in the nervous system obtained general interest among the researchers [12,29]. At the beginning of the twenty-first century, systematic, geometric and mathematical modeling of these phenomena was initiated [7,13,18], and the subject area of chaos synchronization obtained great attention as well [36,40]. Nowadays, the combination of nonlinearity and fractional order dynamics became an interesting research area [2,38].
In general, in biomedical problems the main sources of nonlinearities and variable coupling originate: -From the mass action law due to which the products of various integer or rational powers of concentrations occur in the balance equations, -Nonlinear truncations, because the physical interpretation does not allow negative concentrations, -The limitation of the control signals because a reagent can have only positive ingress rate-it cannot be purely extracted from the stirring tank or the living organism, -From the phenomena of "Input Coupling" [35], meaning that by adding some reagent to the system its other components are inevitably diluted.
Diabetes mellitus (DM) is a chronic disease of the human metabolic system regarding the malfunction in the production and utilization of insulin hormone. Several types of DM exist grouped on the reason of DM [1,11]: -The lack of insulin production classifies Type 1 DM (T1DM); -Resistance against the effect of insulin categorizes Type 2 DM (T2DM); -Both the aforementioned cases create Double DM (DDM); -DM during pregnancy classifies Gestational DM (GDM); -Finally, there are DM caused by genetic disorders.
The most dangerous type of DM is T1DM which occurs when the patients' own immune system identifies the pancreatic β-cells-which producing the insulin hormone-as targets and destroys them during an autoimmune reaction. Because of the lack of internal insulin, the patients need external insulin in order to avoid metabolic collapse. Beside the shortterm handling, maintaining the long-term variability of glycemia is also important to avoid the long-term side effects of the disease [5].
Diabetes is not curable, but treatable. The treatment depends on the type of diabetes; however, in case of T1DM and over time in T2DM this means external insulin administration by insulin pen (manual) or insulin pump (semiautomated) [17]. From engineering point of view, the best solution in order to reach a closeto-normal glycemia during DM treatment is the semiautomated insulin pump therapy, where the electromechanical device administers the required insulin based on a developed control algorithm [4].
Modeling and identifying DM models are not trivial tasks. Almost all available models contain high nonlinearities which make the control designing procedure difficult. The input time signals have impulse nature, since the meal intakes and the external insulin administration can be modeled as impulse functions [11,31]. Moreover, the output of such models can be compared only with quantized blood glucose data, since the commercially available continuous glucose monitoring sensors (CGMS) [4] used with the insulin pump systems measure on every 5 min due to technological limitations.
In the recent years, several advanced control solutions appeared with regard to maintain the glycemia [30]. However, most of them are model-based solutions using different simplifications of the nonlinear problem because of the aforementioned unfavorable circumstances [16,21].
The current research work focuses on a recently appeared robust nonlinear solution, the robust fixed  [32], that does not need exact models just roughly approximation of the real-world problem. The RFPT-based design method is demonstrated on the DM control problem using two control design approaches based on the affine or non-affine model of the physiological problem. The paper is structured as follows. The first section contains the detailed description of the applied diabetes model, while the second section introduces the RFPT-based methodology. In the third section, the controller design procedure is demonstrated followed by the research results. The final section contains the conclusions and future work possibilities.

Diabetes model
In this study, a recently appeared glucose-insulin model is investigated, developed in order to increase the efficiency of identification from real patients' samples [23]. The state-space representation of the model is given as follows:  [23]. The model has two inputs, namely the external insulin infusion rate u(t) (U/h) and the carbohydrate (CHO) intake r (t) (mg/min 2 ), and one output, the glycemia, G(t) (mg/dL) used as a state of the model as well. Other states of the model are the insulinemia, I (t) (U/L) and the digestion of CHO, D(t) (mg/dL/min). The first subsystem (Eq. 1a) is responsible to simulate the glucose dynamics with regard to the external and internal glucose appearance, the effect of insulin and the internal insulin-independent glucose consumption. The second subsystem (Eq. 1b) describes the insulin dynamics including the changing of insulinemia and external insulin intake, while the third subsystem (Eq. 1c) presents the digestion dynamics and creates connection between the CHO r (t) in meal and D(t).
In order to determine the relative order of the necessary control, the order of the time-derivative of G(t) has to be found. This can be immediately set by the control signal u(t), the insulin ingress rate. For this, the "effect chain" of the control signal has to be clarified.

Effect chain of the control signal
According to (1b), u(t) immediately influences thë I (t). SinceÏ (t) occurs in the third time-derivative of G(t), Eq. (1a) has to be differentiated two times: Via substituting (1b) and (1c) into (2b), the control equation can be obtained: ...
2484 L. Kovács from which the necessary control signal u(t) for the prescribed ... G(t) can be calculated: From Eq. (4), it becomes clear that the u(t) control signal directly affects the ... G(t)-that means the control law has to cover this connection.

The robust fixed point transformation-based adaptive control method
The idea of solving nonlinear equations via iterative techniques has long traditions in numerical computing. In a wider context, the original task can be transformed into a fixed point problem that in the next step can be solved by iteration. For instance, the Newton-Raphson algorithm is a classic example that seems to be one of the fundamental methodologies and attracts great attention even nowadays [8,20,28,39]. In the sequel, the transformation of the adaptive control task into a fixed point problem is briefly highlighted. This is followed by the creation and the convergence properties of the iterative control signal.
3.1 Determination of the relative order of the control task: the kinetic tracking error prescription and the "response function" In order to determine the relative order of the control task, the first step is to consider the physical quantity for which a nominal time-variation or nominal trajectory G N (t) is defined in the given task. This step can be made on the basis of purely kinetic considerations by trying to prescribe the appropriate order time-derivative of the controlled quantity that instantaneously can be affected by the control signal. As it will be seen in Sect. 4, in our case, according to the model in use, the third time-derivative of the glucose concentration of the blood ... G (mg/dL min 3 ) can be directly influenced by the control input u (U/h) that is the external insulin infusion rate. For instance, by considering the integrated tracking error defined as , and by introducing a positive real number 0 < (s −1 ), we may wish to have d dt + 4 e int (t) ≡ 0 that yields the "desired system response" as: In the possession of an available approximate system model, in the given situation, the controller can estimate the appropriate value u that, if the model would be exact, just would generate this ...
Due to modeling and/or state estimation errors, when this control signal u is applied on the actually controlled system, the realized (and measurable) "realized response," i.e., ...
On this basis, a "response function" can be defined that for an arbitrary input ... G In yields the realized response as .. .
. . , in which in the place of the symbol "…" the zero-, first-and second-order derivatives of G(t) and the other state variables of the system can be understood. Since ... G can be instantaneously modified by u, while the other arguments in the place of the symbol "…" vary only slowly, we can use the approximation .. .
In the lack of information on the exact model parameters, the analytical expression of f is not available for the controller. However, the pairs made of ... G and ... G In are always known: the input value is determined by the controller, and ... G is measurable. In the sequel, by the use of the response function, an iteration is suggested to find the appropriate value 3.2 Transformation of the control task into a fixed point problem: the "robust fixed point transformation" Assume that we have a digital controller, and in each control step we can make exactly one step of iteration by the use of a function H defined as follows: ...
where the real numbers K c , A c , and B c are the adaptive control parameters. This simple function was introduced in [33]. If ... that cannot be used for control purposes. This construction evidently corresponds to the requirement of causality; the controller learns from the experience made in the "recent past": the signal to be used in control cycle number (n + 1), i.e., ... In the next subsection, the convergence properties of the here defined sequence are considered.

Guaranteeing the convergence of the iteration in general
As is well known, a Banach space (as a set, casually denoted by B) by definition is a complete, linear, normed metric space [14], i.e., it has the following properties: -Linearity: ∀α, β ∈ C and x, y ∈ B, the linear combination is defined and belongs to the space: αx + βy ∈ B; -Existence of a norm for defining a metrics: ∀x ∈ B ∃ x ≥ 0 so that from a = 0 it follows that a = 0 (i.e., the norm separates points), ∀α ∈ C αx = |α| · x (absolute scalability), and x + y ≤ x + y (norm inequality); by the use of this norm, the metrics or the distance between the elements as ρ(x, y) def = x − y can be defined; -Completeness: By the use of the concept of the norm the so-called Cauchy-sequences can be defined as follows: a sequence {x n ; n ∈ N} is a Cauchysequence if ∀L ∈ N x n+L − x n → 0 as n → ∞; completeness means that each Cauchy-sequence must be convergent in a complete space, i.e., for the above sequence ∃x ∈ B so that x n − x → 0 as n → ∞.
By the use of the norm, for the functions Φ : B → B, contractive ones can be defined as follows: By the use of a contractive function, Cauchy-sequences can be generated in the following manner: This sequence is evidently a Cauchysequence since: Due to the completeness of B, ∃x ∈ B so that x n − x → 0 as n → ∞. It is easy to show that x is the Fixed Point of Φ, i.e., Φ(x ) = x . By utilizing the properties of the norm, it can be written that: As a result, the advantage of Banach spaces is the use of the above simple and practical argumentation in the case of quite "abstract" and "complicated" sets [3]. For instance, the quadratically integrable functions of modern quantum mechanics form a Hilbert-Space, that is only a special example of Banach spaces. The above simple considerations allowed the application of the fixed point transformation-based approach for the control of systems of R n → R n , n ∈ N-type response functions and made it possible to further clarify the conditions of convergence in [9].

Convergence conditions for SISO systems
In our case, the ... G ≈ f ... G In response function corresponds to f : R → R, so we can use the Banach space of real numbers with the norm x def = |x|. Since our function is differentiable, we can use a simple integral estimation for guaranteeing contractivity: Since |K c | is very big, in this case the initial element of the iteration will be either in the (−K c , ... G ) interval or it will be greater than ... G ; therefore, the iteration will converge to ... G . An advantage of this control method is that it does not require very precise setting of its adaptive control parameters. The actual setting concerns the speed of convergence, therefore, to some extent the precision of trajectory tracking. In the field of life sciences, this fact is very important from practical point of view, because the method is not based on "exact proofs" for which the necessary conditions rigorously have to be guaranteed.

Controller design
In the followings, the controller design method is demonstrated in the given case by using the aforementioned theorems and methods started with the realization of the affine and approximate system models.

The affine model
It is evident that (4) corresponds to an affine structure as the relationship between u(t) and ... G(t) is concerned. It is reasonable to assume that any abrupt jump in u(t) immediately affects the instant value of ... G(t); hence, the "additive" parts of the affine model vary only slowly. As a result, in the RFPT-based control design the actual value of the control sequence r n will be the "required" third derivative of G. It will be referred to as ... G(t)Req in the sequel. According to the available model to ... G(t)Req, the control signal u(t) Req becomes as follows: The phenomenological restrictions that are so typical in the control of T1DM obtain significance at this point.
Practically only G(t) can be measured by an appropriate continuous glucose monitoring sensor with 5 min cycle time. No direct measurement possibilities exist for measuring D(t) and I (t) in the practice. However, in principle r (t) may be known, as it depends on the action of the patient, but it cannot be expected that the patient "manually" provides the controller with this information. Therefore, it is assumed that the controller can "detect" the CHO intake through observing some increase in G(t).
Consequently, it can be stated that in practice there is no viable way to obtain information on the actual value of the additive parts of the affine model. The main feature of the RFPT-based adaptive controller that can work with an incomplete model and approximate it well fits to this practical problem: we do not need the application of complicated state estimators to estimate this term. In our model, this unknown contribution is denoted as an "AffineAdditive" constant term as follows: In this approach, the information on the variables D(t) and I (t) is completely neglected.

Approximate model
As no real measurements can be done for the estimation of the actual D(t) and I (t) values of the patient, an alternative possibility is the application of some "approximate model" (its parameters are denoted by the symbol ∼) to estimate them by solving the complete equations of motion for the approximate model taking the same control signal as the actual patient u(t), and producing the drift of "approximate state variables"Ġ(t),İ (t),Ḋ(t). In the estimation of u(t) Req the "approximate quantities"Ĩ (t),İ (t),D(t),Ḋ(t) are substituted from a computer program that emulates the behavior of the approximate model, but it takes the measurable actual G(t) value and-in the lack of information-instead of the actual input r (t) it takes zero: This approach may give less work to the adaptive compensation than that of the simple "affine model."

Control law
During the study, a kinematic-type PID-control law was used. This is an appropriate choice if the goal of the control is trajectory tracking as it is in this case. Since the control signal affects the third derivative of the variable to be regulated, the control law of the same order has to be used: which determines the following desired ... G Desired (t) function: where the G N (t) is the reference (nominal) blood glucose (BG) level, G(t) is the actual (real) BG level, and the error is the G N (t) − G(t) that has to converge to zero over time.

Final control environments
The final control environment consists of the PID-type kinematic prescriptions (the control law), the adaptive block and the affine model (Fig. 1). However, in this study another realization possibility was investigated as well, where an approximate parallel simulated model provides the non-measurable estimated states. The structure of it is presented in Fig. 2. Note that the PID-type prescription and the adaptive block were the same.

Results
In order to test the controller in the case of unfavorable circumstances, long-term simulations were applied and   Fig. 3 Used carbohydrate CHO intake a unique glucose input function was designed. In every case, the goal was that the controller should provide an appropriate control signal by which the glycemia of a patient can be stable and appropriate in long-term beside continuous glucose disturbances.

CHO intake
The glucose input specificity of the used model required a special input function. The designed function consists of an additive mixture of an arbitrarily selected sinusoidal disturbance signal and an impulse kind signal, r (t): where the s h is the impulse height and s w is the impulse width of the signal. The model needs the derivative of the designed function dr (t)/dt as CHO input. Figure 3 shows the output of the designed specific CHO input function.
The primary goal of this study was to prove the usability of the RFPT-based controller design opportunity in case of the T1DM model created in order to ease the identification procedures. The method requests approximated models instead of exact patient models that allows using one of the parameter sets from the study of the model belonged to an identified patient (Patient 1, [23]) given in Table 2.
The AffineAdditive elements were set to zero during the simulation. Naturally, the parameters from Table 2  can only be used in the affine case. In the approximate case, the values of the state variables are nonmeasurable; however, roughly estimable from parallel simulation enough in order to efficiently use the RFPTbased method. The parameters in this case were half of the original values; namely, except the body weight, every parameter in the approximate case was equal with 0.5 times of its original value. Beside the selection of the used model parameters, the appropriate selection of the control parameters is also important. The general RFPT-based controller parameters are those connected to the adaptivity, namely A c , B c and K c . Moreover, depending on the applied control law, different further control parameters may occur. In this study, a kinematic control law was used (see Sect. 4.3) with one tunable gain parameter. The values of the adaptivity parameters were adjusted to the magnitudes of the controlled variable (the third derivative of G(t)). The last selectable variable is the reference (nominal) BG level G N which is used in the adaptivity block and the control law as well. Due to the desired goal, to prove the usability of the RFPT-based controller, the same control parameters and reference BG level were used during the simulations, without online parameter tuning. However, in other applications these properties of the RFPT-based controller design were successfully tested [34,37]. The selected control variables of the T1DM case is given in Table 3.
The arbitrarily selected simulation length was 4000 min (more then 66 h, or almost 3 days), which is enough to demonstrate the benefit of the RFPT-based controller, i.e., the controller adapts to the patient's needs.  Figure 4 presents the blood glucose level over time in the affine case. It should be noted that the desired blood glucose interval is 70-120 mg/dL (the healthy human blood glucose interval). Over 120 mg/dL hyperglycemia (high BG), under 70 mg/dL hypoglycemia (low BG) is diagnosed. The latter is the most dangerous for a T1DM patient and should be completely avoided. Hyperglycemia, however, could be tolerated, but the amplitude should be reduced to 140 mg/dL.
Evaluating the results presented in Fig. 4, it can be seen that after the first transient the controller reacts to the increasing BG level and administers insulin in order to avoid hypergycemia. As a result, the BG level finally reaches the nominal BG level G N ; however, the controller was continuously operating to prevent the unfavorable effects.
In the non-affine case of Fig. 5, the controller acts faster, since the simulated approximate model signals are available and the controller has direct, but roughly approximated information about the possible internal states. The mild waviness in the figures is the effect of the control (insulin) signal, while the peaks occurred are discussed later. Figures 6 and 7 show the tracking errors in the affine and non-affine cases. The same conclusions can be observed; namely, in the affine case the tracking error decay is slower (as the RFPT-based controller did not have only indirect information), but the controller works efficiently over time. The waviness effect appears here as well, since the error signal stands from G N − G realized (t). Figure 8 shows the injected insulin over time in the affine case. Due to the affine model's structure (Eq. 9) insulin peaks occur over time because the third derivative of the required G(t). Originally, these effects come from the food intake signals and reflects in the ... G(t)Req signal, respectively.  The same effects can be seen in the non-affine case (Fig. 9). The average magnitude of the peaks is almost the same. The main difference is that the controller has approximated internal information about the states. This knowledge insulinemia that allows us to have a nonzero initial control signal and have the non-smooth insulin signal around the peaks. In Figs. 10 and 11, one can see the output of the realized (the real system's) answer, desired (based on the control law) and required (recommendation of the adaptivity block) third derivatives of the BG level G(t) in both the affine and non-affine cases. In Fig. 10, the realized signal tends to the desired signal after the diversion caused by the insulin signals; hence, the adaptation works well. In Fig. 11 the desired and required signals are almost the same, which is the direct con-  Fig. 9 Injected insulin in the non-affine case sequence of that fact that the controller has internal, however, roughly approximated information about the states. The realized ... G(t) was tended to the desired signal as the low magnitudes and due to the fact that the desired and required signals were almost the same. Figure 12 shows insulinemia variation I (t) over time. It can be seen that the control insulin signal results a stable internal insulin level, nonetheless the presence of the insulin peaks. The delays of the effect coming from the model's structure are also visible, since despite  Fig. 11 Desired, required, realized third derivatives of G(t) in the non-affine case the immediate insulin signal from the beginning the insulinemia initially decays, but it stabilizes over time.

Injected insulin
In the non-affine case (Fig. 13), the waviness appears come from the insulin peaks.

Conclusion
In this paper, the usability of RFPT-based controller design method was reported in the case of a type 1 diabetes nonlinear model optimized to long-term identification purposes, but used for control purposes as well.
In line with the requirements of the T1DM model, an appropriate feed intake function was designed and long simulation time was used in order to provide unfavorable circumstances to test the behavior of the developed RFPT-based controller. During this study, fixed control parameters were used without advanced optimization techniques. On these considerations, the applicability of the novel RFPT method has been demonstrated. With online optimization techniques, better performance is expected that is a next step of the presented research.