HNLN: Artificial Pancreas
The hybrid automaton of an artificial pancreas is presented as below. The system has 10 variables whose evolutions are governed by the following ODE.
X' = - 0.0278 * X +0.0278 * ( 18.2129 * Ip - 100.25)
Isc1' = - 0.0171 * Isc1 + 0.97751710655*InsulinRate
Isc2' = 0.0152 * Isc1 - 0.0078 * Isc2
Gp' = 4.7314 - 0.0047 * Gp - 0.0121 * Id - 1 - 0.0581*Gp + 0.0871*Gt + Meal
Gt' = -0.0039* (3.2267+0.0313*X) * Gt * ( 1 - 0.0026 * Gt + (2.5097e-6) *Gt^2) + 0.0581 *Gp - 0.0871* Gt
Gs' = 0.1*( 0.5221 * Gp - Gs)
Il' = - 0.4219 * Il + 0.2250* Ip
Ip' = - 0.3150 * Ip + 0.1545 * Il + 0.0019 * Isc1 + 0.0078 * Isc2
I1' = -0.0046 * ( I1 - 18.2129 * Ip)
Id' = -0.0046 * (Id - I1 )
The input Meal is defined by a peice-wise polynomial function:
0 <= t < 30: (1.140850553428184e-4)*t^2*uc + (6.134609247812877e-5)*t*uc
30 <= t < 80: (5.252705445621361e-5)*t^2*uc - 0.007468661585336*t*uc + 0.281302431132355*uc
80 <= t < 360: (1.245404770117318e-7)*t^2*uc - (9.112517786233733e-5)*t*uc + 0.026475608001390*uc
360 <= t < 400: -(6.307316106980570e-5)*t^2*uc + 0.048262107570417*t*uc - 9.190266060911723*uc
400 <= t < 500: (3.552954405332104e-6)*t^2*uc - 0.003423642950746*t*uc + 0.823855671531609*uc
t >= 500: (1.113158675054563e-8)*t^2*uc - (1.482047571340105e-5)*t*uc + 0.004900138660598*uc
such that uc is a time-invariant uncertainty ranging in [70,80].
The insulin rate is updated according to the following rule:
1) If the glucose level Gs is increasing and reaches 80, then we set insulinRate = 0.1.
2) If the glucose level Gs is increasing and reaches 120, then we set insulinRate = 0.2.
3) If the glucose level Gs is increasing and reaches 180, then we set insulinRate = 0.5.
4) If the glucose level Gs is increasing and reaches 300, then we set insulinRate = 1.4.
5) If the glucose level Gs is decreasing and reaches 75, then we set insulinRate = 0.05.
6) If the glucose level Gs is decreasing and reaches 115, then we set insulinRate = 0.1.
7) If the glucose level Gs is decreasing and reaches 175, then we set insulinRate = 0.2.
8) If the glucose level Gs is decreasing and reaches 295, then we set insulinRate = 0.5.
We consider the time horizon [0,720] and jump depth 8. The initial set is given by
X in [0,0]
Isc1 in [72.43, 72.43]
Isc2 in [141.15, 141.15]
Gp in [229.824, 268.128]
Gt in [162.45, 162.45]
Gs in [120, 140]
Il in [3.2, 3.2]
Ip in [5.5, 5.5]
I1 in [100.25, 100.25]
Id in [100.25, 100.25]
The unsafe set could be Gs >= 300. Also we would like to ensure that the glucose level should be stabilized between 240 and 280 over the time interval [600,720].
We created a SpaceEx model for the artificial pancreas system descibed in
http://www.cs.colorado.edu/~srirams/papers/arch-2017-multi-basal.pdf, which is similar but slightly different to the artificial parncreas system above. The SpaceEx model as well as the results from reachabilty analysis with CORA are attached to this post. For reachability analysis we used the following initial set:
X in [0,0]
Isc1 in [72.43, 72.43]
Isc2 in [141.15, 141.15]
Gt in [162,45, 162.45]
Gp in [229.824, 268.128]
Gs in [120, 140]
Il in [3.2, 3.2]
Ip in [5.5, 5.5]
I1 in [100.25, 100.25]
Id in [100.25, 100.25]
Gs in [120, 140]
t in [0, 0]
Further we used uc in [70,80] for the uncertain system input. We calculated the reachable set for the time interval t in [0,720] min.
We wanted to share the results we obtain with CORA for this benchmark. The code is compatible with the current CORA version on the public repository https://github.com/TUMcps/CORA
There seems to be an inconsistency between the model used in the paper and the model we use, since for our case the simulation violates the specification G_s < 300.
