Receding Horizon Control of Type 1 Diabetes Mellitus by Using Nonlinear Programming

1Doctoral School of Applied Informatics and Applied Mathematics, Óbuda University, Bécsi Street 96/B, Budapest 1034, Hungary 2Mathematical Sciences Research Center, Karachi, Pakistan 3Antal Bejczy Center for Intelligent Robotics (ABC iRob), Óbuda University, Bécsi Street 96/B, Budapest 1034, Hungary 4Physiological Controls Research Center, Óbuda University, Bécsi Street 96/B, Budapest 1034, Hungary


Introduction
The advanced control solutions have inevitable role in today's medical practice regarding the control of physiological processes [1].Many control solutions are under development that can be used for various kinds of control problems.Advanced control methods have been successfully applied for physiological regulation problems, for example, control of anesthesia [2,3], angiogenic inhibition of cancer [4,5], immune response in presence of human immunodeficiency virus [6], and regulation of blood glucose (BG) level [7][8][9][10] as well.
Diabetes Mellitus (DM) is the collective name of several chronic diseases connected to the metabolic system of the human body.In most of the cases, the DM condition appears due to the issues related to the insulin hormone [11].The insulin is the key hormone which makes the glucose molecules possible to enter from the blood into the glucose consuming cells through the insulin-dependent gates on the cell-wall [12].
There are many types of DM.The most dangerous is the Type 1 DM (T1DM) where the metabolic system is not able to function normally due to the lack of insulin.Type 2 DM (T2DM) is the most widespread kind of DM and it occurs mostly because of the lifestyle.In this case usually the blood glucose and insulin levels continuously increase over a long period of time.Due to the extreme glucose and insulin load the cells become resistant to the insulin over time.In order to compensate this condition the body produces more and more insulin that leads to the "burnout" of the pancreatic -cells that produce the hormone.At this point the T2DM turns into T1DM.Other frequently occurring type is the Gestational DM (GDM) from which women may suffer during pregnancy.Usually, this condition is 2 Complexity temporary; however, sometimes it turns into T2DM and becomes permanent [13][14][15].
In case of DM the application of these kinds of advanced control techniques has high importance.Due to the nature of the phenomenon to be controlled the researchers on the field have to face many challenges such as high nonlinearities, model and parameter uncertainties, and even time-delay effects, as well.However, regardless the type of DM a few common control goals can be defined: keeping the glycemia (the BG level) in a the healthy range; totally avoiding the hypoglycemic periods; and avoiding the high BG variability as much as possible [16][17][18].
In this research we have investigated the T1DM.As we mentioned, T1DM is the most dangerous condition because the patients need external insulin intake in order to keep their metabolic status on appropriate level.In this case the patient's pancreatic -cells are terminated by the immune system of the patient during an autoimmune reaction.As a consequence, these patients are not able to produce insulin-which leads to short term starving, coma, or even death [11].Furthermore, the appropriate treatment-how the insulin is administered-is also important to avoid long term side effects, for example, the chronic failure of peripheric vasculature [13].
By using advanced control techniques not just acceptable control action but higher treatment quality can be obtained.The selected control methods have to handle the already mentioned unfavorable effects-such as the nonlinearities and so on-as well.In case of T1DM many solutions are available; however, all of them have their own limitations, simplifications, and restrictions-thus, none of them are general [10].In these days from control point of view the most beneficial approach is the Artificial Pancreas (AP) concept.This idea aims to imitate the regular operation of the pancreas from the insulin production point of view, namely, administering insulin demands on the needs determined by the BG level [19].Thus, we have to face contradictory requirements: the generalization and personalization as well.
One of the mostly used algorithms is the modified proportional-integral-derivative (PID) solutions due to their simplicity and flexibility.Moreover, several clinical trials have been done by using this methodology and investigate its effectivity [20][21][22].Linear Parameter Varying (LPV) modelbased solutions have high importance, since the uncertainties can be handled with high efficiency by them [7,9,23].The Tensor Product (TP) model-based techniques also represent interesting directions, since they can be combined by Linear Matrix Inequality (LMI) based control and LPV methodology as well [24,25].
The most frequently used method is the Model Predictive Control (MPC) regarding the control of DM [19,[26][27][28].MPC is a widely used approach in various research fields as well [29][30][31][32].In general the goals of the MPC applications are tracking and stabilization [33].In case of an MPC the actual value of the control signal () is obtained by solving an openloop control problem over a finite horizon.Of course, some feedback is present because the starting point of the next horizon is the last realized point of the previous one.The type of the optimal control problem depends on the type of the control task.The classical MPC is realized in a framework of cost-function-based optimal control where the dynamics of the system to be controlled can be considered as a set of constraints.Furthermore, the cost function regularly contains terms which depend on the tracking error and the control signal itself.Optionally, a cost contribution that depends on the tracking error at the end of the horizon can be added as well.
In case of the Receding Horizon Control (RHC) via Nonlinear Programming (NP) based approximation the state variables and the control inputs of the system are considered over a discrete time grid.At each point of the grid the Lagrangian multipliers determine the reduced gradient which is driven into zero numerically to find the optimal solution.This solution consists of the estimated values of the state variables and control signals over the finite horizon.In that case if the dynamic model of the system is not precise, the optimal design can be used only for consecutive finite horizons since the actual state of the controlled system propagates according to its exact dynamics.In order to minimize the effects coming from the inaccuracies the actual measured state variable from the end of the previous horizon is applied as starting point for the next one in the next horizon-length design [34].The RHC framework can be hardly combined with the Lyapunov function based control.However, certain approaches can be found in the literature where the Lyapunov stability [35] and RHC were successfully combined for specific cases [36][37][38][39].
Alternative solutions also exist which can be used instead of Lyapunov's stability theorem.The Robust Fixed Point Transformation (RFPT) based control [40,41] uses Banach's fixed point theorem [42] to transform the control problem into a fixed point problem which can be solved iteratively.This method allows designing a robust iterative adaptive controller which can avoid the main limitations of RHC if these are combined.
In this study we present the first part of our research, namely, the design of an appropriate RHC controller on NP basis which can be completed by RFPT in our further work.
The paper is structured as follows.First, the applied diabetes model is introduced.After that, the RHC design based on the NP approach is presented.Then, the results are introduced with their discussion.Finally, we conclude our work and provide an outline regarding our further research.

Type 1 Diabetes Mellitus Model
During our research we have applied a modified Minimal Model [43] which originates from the model of Bergman [46].This model has several beneficial properties, such as simplicity, good transformability, and flexibility and it is based on simpler biological considerations.The main goal of the model is to describe the glucose-insulin dynamics, namely, to define the connection between the blood glucose and insulin levels.However, in order to characterize the daily life of a T1DM patient this model has to be extended with additional submodels.These submodels are the absorption of the external glucose and insulin intakes.During the daily routine these substances are not directly injected to the blood stream-however, this can occur in case of persistent hospitalization.Instead, the carbohydrate is consumed via food intake and the insulin is entered through the extracellular tissue matrix under the skin [13].Thus, their appearance does not contain sharp peaks: it happens through longer dynamics.
The glucose and insulin absorption are described by ( 1)-( 4), respectively.These submodels originate from the Cambridge model [44], but we applied them in appropriate dimensions to insert them into the core model.The core model is described by ( 5)- (7).
The state variables in (1)-( 7) have the meaning and purpose as follows. 1 () mg/dL and  2 () mg/dL are the primary and secondary compartments belonging to glucose, where the time constant   determines how long it takes for the meal to be absorbed after consumption in time. 1 () mU/L and  2 () mU/L are the primary and secondary compartments belonging to insulin, where the time constant   determines how long it takes the insulin to be absorbed after injection (to the extracellular space) in time.Variable () mg/dL is the blood glucose (BG) concentration-the so-called glycemia-and () mU/L is the blood insulin concentration and () 1/min is the insulin-excitable tissue glucose uptake activity-which describes the connection between the blood's glucose and insulin levels, respectively.
From system engineering point of view the external glucose, namely, the food intake, can be handled as disturbance.In this case () g/min is the disturbance input.It can be inserted to  1 () via ((1000  )/(    )) complex which describes the bioavailability of the glucose from complex carbohydrates.The control signal () mU/min-the injected insulin-is directly connected to  1 ().More detailed description of the used model parameters can be found in Table 1 and in [43][44][45].

A Nonlinear Programming Approach with regard to the Receding Horizon Controller
3.1.Nonlinear Programming in General.The numerical approximation of a given problem can be considered in the following way: (I) Determination of a discrete time grid which is dense enough from the given application point of view as { 0 ,  1 =  0 + Δ, . . .,  +1 =   + Δ, . . .,   }, where  ∈ N.Here  0 belongs to the initial and   belongs to the final time instant of the motion to be taken into account, respectively.
(II) Assume that the nonlinear equation of the system to be controller is ẋ () = (x(), u()), where x() ∈ R  is the state and u() ∈ R  is the input vector, respectively.Furthermore, the x 0 ≡ x( 0 ) initial condition of the system is given.
(III) Consider that x  (  ) is the nominal trajectory to be tracked by the states of the system over time.Here x  (  ) ≡ x   in the given time grid determined by (I).(IV) The nominal prescribed trajectory cannot be exactly realized during the execution of the control task.However, several restrictions-in order to enforce the realization of this trajectory-can be prescribed by using a predefined (x(), u()) cost function in each point of the grid.
(V) (x(), u()) ≥ 0 is able to express several requirements, although these are often contradictory.It can be constructed as the sum of nonnegative terms which can be the differentiable functions of the control signal and the state variables-expediently.The drastic control signals can be avoided as well by prohibiting them via the cost function.
(VI) Terminal conditions can be embedded into the prescription depending only on x  in the last step   .
(VII) During an optimal control approach, ∑ −1 =0 (x  , u  ) + Φ(x  ) has to be minimized, where Φ(x  ) is an extra weight that belongs to the last point of the trajectory.
(VIII) The cost function ∑ −1 =0 (x  , u  ) + Φ(x  ) cannot be arbitrarily minimized due to the specificities of the system to be controlled.In the minimization the state propagation equation must be considered as a constraint.In order to process this constraint the Lagrangian multipliers can be used as follows.
(IX) ẋ () can be expressed from the state propagation equation and as a numerical approximation (x +1 − x  )/Δ ≈ (x  , u  ).By using this estimation, the control task can be expressed in the following way: and { 0 , . . .,  −1 } are the Lagrangian multiplierswhich are used in accordance with the optimization task to be solved by the reduced gradients method.

The Applied Cost Function.
During the development of the appropriate cost function-which fits to the given problem-the specificities of model ( 1)-( 7) should be taken into account.
The main limitation coming from the model is the amount of injectable insulin and the fact that we have only one control signal.The control signals at the grid points of the horizon are independent variables in the optimization problem.However, the application of a specific form of limitation on them-a "bias"-is reasonable.Thus, the control signal in this construction should be limited in accordance with the phenomenon to be controlled.Instead of the control signal itself, another variable should be selected as independent variable to avoid the initial value problems causing the rough numerical approximation at the beginning of the optimization.This is caused by the high nonlinearity in the model.
Another property of the model is that only the blood glucose level () can be measured.Thus, only this state variable can be embedded into the cost function to be developed.We do not have internal information about other state variables of the process to be controlled.
Accordingly we have applied a more specific form of the (8) cost function: where   =  bias + tanh(V  ).
We have developed a strongly nonlinear cost function in which all requirements can be embedded against the control action to be reached during control.

𝐽 (𝐺, 𝑢)
The tracking error in ( 10) and ( 11), namely, the deviation of the realized blood glucose level () from the nominal blood glucose level   (), can be calculated as    −   at the grid points.The absolute value of this difference can be determined by parameters  1 and  3 that contribute the belonging level of "penalty" prevailing in the cost function.In that case if  1 > 1 and  3 > 1 beside |  − | <  1 and |  final −  final | <  3 , the contribution to the cost function is low.However, if |  − | >  1 and |  final −  final | >  3 then due to the power terms the contributions to the cost of these terms are drastically increasing.The  1 and  3 weighting terms can be used for different purposes. 1 =  3 = 1 provides proportional contribution; namely, the pure deviation will be better prohibited.However, if  1 ,  3 < 1, then the smaller deviations will be prohibited relatively better than the bigger ones.The role of  2 and  2 parameters are similar; namely, the applied control input can be prohibited by applying them.The  parameter allows modifying the enforcement of the effect of the control signal in the cost function.

Applied Nonlinear Programming Solver.
In order to implement the reduced gradients based optimization as Nonlinear Programming task various software products can be used.Due to its complex embedded nonlinear solvers and the easy-to-use property we have selected the Microsoft's EXCEL (MS EXCEL).The MS EXCEL's "Solver" module is produced by an external firm Frontline Systems, Inc.There are various solutions implemented into the Solver module, including the Generalized Reduced Gradient (GRG) method which can be used in the given case as well.The GRG is based on [47,48] and its usability has been proved in various fields of research.For example, it was successfully applied for macroeconomic optimization problems [49], operational research problems related to economics [50], optimization problems regarding transportation [31], decision prediction models [51], optimal control problems in continuous time domain [52], and so on.The optimization framework was built up by using the "Visual Basic" packages of the software "MS EXCEL 2016" under "Windows 10" operating system.The functional dependencies have been constructed in Visual Basic, and the grid points, parameters, and so forth have been defined on dedicated Worksheets-by using these data the Solver can be easily set up.

Results
In order to test the realized control framework we investigated two scenarios.In the first scenario 25 g CHO was considered-5 g over 5 minutes in each 240th time instant from the 60th one.
In the second scenario we considered 50 g CHO intake 10 g over 5 minutes in each 240th time instant from the 60th time instant.In both cases 200 time horizons have been considered within 10 grid points; thus the total simulated time domain was 2000 minutes.The resolution Δ was 1 minute in accordance with the properties of the model.
The applied cost-function parameters (which represent the control parameters in this regard) can be found in Table 2.
It should be noted that we have used permanent reference trajectories in both cases denoted by   = 90mg/dL in ( 9)- (11).
Therefore, in accordance with the aforementioned details, the goal of the control becomes to keep   = 90 beside respecting the predefined  bias .In this manner-via the cost function-the deviations from these predefined values have been "punished"; namely, the value of the cost function became higher.
During the examinations we have considered the following range of blood glucose as "healthy" in accordance with [11]: from 70 to 180 mg/dL.

Results of Scenario 1.
In the following the results of Scenario 1 are presented.First, the disturbance signal is shown by Figure 1.The applied-calculated-control signal can be seen in Figure 2. It is clear that the controller is able to administer the insulin in accordance with ( 9)- (11) where   =  bias + tanh(V  ) is prevailed in the control action.
Figures 3 and 4 show the absorption of glucose and its appearance in the blood with the dynamics determined by the model.Figures 5 and 6 show the absorption of insulin from the interstitium and its appearance in blood.The main point can be seen on Figure 7.The controller is able to satisfy the determined conditions and the BG level (()) is inside the predefined range-no hypo-and hyperglycemia occurred.The BG level approaches the selected reference trajectory (  ) as it is expected.
On Figures 8 and 9 the variation of blood insulin level and the intermediate state variable can be seen.() determines how the blood insulin level affects the blood glucose level, namely, the connection between them.

Results of Scenario 2.
The applied disturbance input in accordance with the detailed protocol can be seen in Figure 10.In this case, we have applied higher inputs in order to be sure that the developed control framework is able to deal with unfavorable external excitation.
Figure 11 reveals the calculated and administered control signal.As it can be seen, its dynamics are significantly different from that of the previous case due to the different settings in the applied cost function.
Figures 12 and 13 represent the absorption of the glucose and its appearance in the blood with the dynamics determined by the model.Figures 14 and 15 show the absorption of insulin from the interstitium to the blood.
The main result can be seen in Figure 16.Though we have drastically increased the disturbance input signal the controller was able to deal with the situation and realized appropriate control action-without domain violation from the determined  bias point of view.The blood glucose level is inside the selected healthy range without any hypo-or hyperglycemia.Moreover, the BG level oscillated around the reference trajectory-  -as it was expected.
Figures 17 and 18 represent the variation of the blood insulin level and the intermediate state ().Due to the higher frequency of the control signal these are oscillating with a higher frequency as well-which is directly reflected in the blood glucose level as well, since the () mediates the insulin's effect on ().

Discussion
In this research our main goal was to design a RHC which is able to control the given patient model and, further, to do that by developing such a RHC controller which can be effectively combined with the RFPT principles in our later work by using the detailed NP environment.These goals have been reached.In our first scenario we aimed to test the control environment under strict constraints-the collection of them can be found in Table 2.In detail,  1 ,  3 = 5 mg/dL and  1 ,  3 = 4 mean that we punished the deviation of () from   heavily above 5 mg/dL during the control action and in the terminal points of the cycles.We applied not just a restriction on the value of () via  bias , but also we penalized the control signals above 2 mU/min as it is reflected in  2 and  2 .
In the second scenario we investigated the "robustness" of the controller by using highly oscillating, unfavorable disturbance input in a cyclic way.As in the previous case, the controller was able to perform well-beside all of the restrictions that have been considered.
The main results can be seen in Figures 7 and 16.It is clearly visible that the main requirement has been satisfied, since the BG level was kept by the controller in the healthy range.

Conclusions and Future Work
In this paper we have proposed a possible design scheme for RHC controller by using a NP approach in order to control Type 1 Diabetes Mellitus.We have applied the MS EXCEL's embedded Solver solution which is based on the Generalized Reduced Gradient method in order to realize the control environment.Two different scenarios have been applied to test our approach-the results have been satisfactory.
In the first scenario we applied "soft" disturbance and smaller penalties via the developed cost function in order to make sure that the controller design is possible and appropriate control action can be achieved by using the continuous optimization.In the second test scenario we used unfavorable, cycling disturbance signal with high amplitude to test the "robustness" of the proposed controller.The developed RHC controller was able to handle the load and provided satisfactory control action.Furthermore, in both cases the BG level was kept in the predefined healthy range.
In our future work we are going to extend the developed RHC controller with an additional RFPT framework in order  to empower it with adaptive property.We will investigate the solution from robustness, adaptivity, and other aspects' points of view.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Calculated control signals.(a) represents the whole time horizon.(b) shows a piece of the whole time horizon between 0 and 40 minutes.

Figure 5 :Figure 6 :
Figure 5: Variation of the first state of the insulin absorption subsystem.

Figure 11 :Figure 12 :Figure 13 :
Figure 11: Calculated control signals.(a) represents the whole time horizon.(b) shows a piece of the whole time horizon between 0 and 40 minutes.

Figure 14 :Figure 15 :
Figure 14: Variation of the first state of the insulin absorption subsystem.

Figure 16 :Figure 17 :Figure 18 :
Figure 16: Variation of the blood glucose level over time.