dreal / dreal2

Please check dreal4 instead.
https://github.com/dreal/dreal4
GNU General Public License v3.0
13 stars 15 forks source link

UNSAT for model with non-deterministic parameter but SAT for deterministic model #57

Open shmarovfedor opened 10 years ago

shmarovfedor commented 10 years ago

I am working with the following model:

#define w 100 //[kg] weight of the patient
#define ka1 0.006 //[min^(-1)] deactivation rate constant (Hovorka 2004)
#define ka2 0.06 //[min^(-1)] deactivation rate constant (Hovorka 2004)
#define ka3 0.03 //[min^(-1)] deactivation rate constant (Hovorka 2004)
#define kb1 0.0034 //[min^(-2) / (mU / L)] activation rate constant (Sriram)
#define kb2 0.056 //[min^(-2) / (mU / L)] activation rate constant (Sriram)
#define kb3 0.024 //[min^(-2) / (mU / L)] activation rate constant (Sriram)

#define tmaxI 55 //[min] time-to-maximum insulin absorption (Hovorka 2004)
#define VI (0.12 * w)//[L * kg^(-1)] insulin distribution volume (Hovorka 2004)
#define ke 0.138//[min^(-1)] insulin elimination from plasma (Hovorka 2004)

#define k12 0.066 //[min^(-1)] transfer rate constant from the non-accessible to the accessible compartment (Hovorka 2004)
#define VG (0.16 * w)//[L * kg^(-1)] distribution volume of the accessible compartment (Hovorka 2004)
#define G (Q1 / VG) //[mmol/L] glucose concentration (Hovorka 2004)
#define F01 (0.0097 * w) //[mmol * kg^(-1) * min^(-1)] non-insulin-dependent glucose flux (Hovorka 2004)
#define Fc01 (F01 * G/(G + 1) / 0.85) //total non-insulin-dependent glucose flux corrected for the ambient glucose concentration (Hovorka 2004)
//#define FR (0.003 * (G - 9) * VG) //renal glucose clearance above the glucose threshold of 9 mmol/L^(-1) (Hovorka 2004)
#define FR 0
#define EGP0 (0.0161 * w) //[mmol * kg^(-1) * min^(-1)] endogenous glucose production 
#define UG1 8 //glucose absorption rate (Hovorka 2004)
#define DG 20 //amount of carbohydrates digested (Hovorka 2004)
#define AG 0.8 //carbohydrate bioavailability (Hovorka 2004)
#define tmaxG 40 //[min] time-of-maximum appearance rate of glucose in the accessible glucose compartment (Hovorka 2004)

[-1.0, 1.0] u;
[-10, 10] x1; //[min^(-1)] (remote) effect of insulin on glucose distribution (Hovorka 2004)
[-10, 10] x2; //[min^(-1)] (remote) effect of insulin on glucose disposal (Hovorka 2004)
[-10, 10] x3; //[min^(-1)] (remote) effect of insulin on endogenous glucose production (Hovorka 2004)
[-10, 10] I; //[mU/L] plasma inslulin concentration (Hovorka 2002)
[-50, 50] S1; //absorption of subcutaneously administered short-acting (e.g. Lispro) insulin (Hovorka 2004)
[-50, 50] S2; //absorption of subcutaneously administered short-acting (e.g. Lispro) insulin (Hovorka 2004)
[-1000, 1000] Q1; //[mmol] mass of glucose in the accessible compartments (Hovorka 2004)
[-1000, 1000] Q2; //[mmol] mass of glucose in the non-accessible compartments (Hovorka 2004)

[0, 500] tau;
[0, 1000] time;

{
mode 1;

invt:
    //(tau <= 500);

flow:

    //glucose-subsystem
    d/dt[Q1] =  - Fc01 - x1 * Q1 + k12 * Q2 - FR + EGP0 * (1 - x3) + (UG1 * 180) / 1000;
    d/dt[Q2] = x1 * Q1 - (k12 + x2) * Q2;
    //insulin-subsystem
    d/dt[S1] = u - S1 / tmaxI;
    d/dt[S2] = (S1 - S2) / tmaxI;
    d/dt[I] = S2 / (tmaxI * VI) - ke * I;
    //insulin-action-subsystem
    d/dt[x1] = - ka1 * x1 + kb1 * I;
    d/dt[x2] = - ka2 * x2 + kb2 * I;
    d/dt[x3] = - ka3 * x3 + kb3 * I;
    d/dt[tau] = 0.5;
    d/dt[u] = 0.0;

jump:
    (G = 10) ==> @2 (and    
                        (I' = I) 
                        (x1' = x1) 
                        (x2' = x2) 
                        (x3' = x3) 
                        (tau' = 0) 
                        (S1' = S1) 
                        (S2' = S2)
                        (Q1' = Q1)
                        (Q2' = Q2)
                        (u' >= 0.35)
                        (u' <= 0.37)
                    );
}

{
mode 2;

invt:
    //(tau <= 500);

flow:

    //glucose-subsystem
    d/dt[Q1] =  - Fc01 - x1 * Q1 + k12 * Q2 - FR + EGP0 * (1 - x3) + (UG1 * 180) / 1000;
    d/dt[Q2] = x1 * Q1 - (k12 + x2) * Q2;
    //insulin-subsystem
    d/dt[S1] = u - S1 / tmaxI;
    d/dt[S2] = (S1 - S2) / tmaxI;
    d/dt[I] = S2 / (tmaxI * VI) - ke * I;
    //insulin-action-subsystem
    d/dt[x1] = - ka1 * x1 + kb1 * I;
    d/dt[x2] = - ka2 * x2 + kb2 * I;
    d/dt[x3] = - ka3 * x3 + kb3 * I;
    d/dt[tau] = 0.5;
    d/dt[u] = 0.0;

jump:

}

init:
@1 (and 
        (I = 0.03)
        (x1 = 0.03)
        (x2 = 0.045) 
        (x3 = 0.04) 
        (tau = 0.0) 
        (S1 = 4.2) 
        (S2 = 4.0)
        (Q1 = 64.0)
        (Q2 = 40.0)
        (u = 0.0)
    );

goal:
@2 (and (tau = 60) (G <= 10));

There is a non-deterministic variable u defined on the interval [0.35, 0.37] and dReach returns UNSAT. But when I change its value to deterministic one from this interval dReach returns SAT

...
jump:
    (G = 10) ==> @2 (and    
                        (I' = I) 
                        (x1' = x1) 
                        (x2' = x2) 
                        (x3' = x3) 
                        (tau' = 0) 
                        (S1' = S1) 
                        (S2' = S2)
                        (Q1' = Q1)
                        (Q2' = Q2)
                        (u' = 0.36)
                    );
...

But satisfiability in deterministic case should imply satisfiability in non-deterministic case.