On the Verification of Weighted Kripke Structures Under Uncertainty

Authors: Giovanni Bacci, Mikkel Hansen and, Kim G. Larsen

Affiliation: Department of Computer Science, Aalborg University, Denmark
e-mail: {giovbacci,mhan,kgl}@cs.aau.dk

We considered the problem of checking weighted CTL properties for weighted Kripke structures in presence of imprecise weights. We consider two extensions of the notion of weighted Kripke structures, namely (i) parametric weighted Kripke structures, having transitions weights modelled as affine maps over a set of parameters and, (ii) weight-uncertain Kripke structures, having transition labelled by real-valued random variables as opposed to precise real valued weights.

Utilities

Quantifier elimination for the theory of linear arithmetic

In this tutorial we will use quantifier elimination for checking the robustness of a WKS when uncertainties on its weights are quantified by means of an absolute error. Mathematica implements a built-in quantifier elimination algorithm, however if one needs a more efficient method one could use Mjollnir.

Mjollnir is a tool developed by David Monniaux that implements quantifier elimination algorithm for the theory of linear real arithmetic. The tool is freely available at http://www-verimag.imag.fr/~monniaux/download/Mjollnir_2012-02-24.zip

Once installed you shall set the installation directory as done below.

Options[QE] = {Path-> "/Users/giovbacci/AAU/QuantifierElimination/Mjollnir/trunk/",Method->"Mjollnir"};
QE[expr_,OptionsPattern[]]:=
If[MatchQ[OptionValue[Method],"Mjollnir"],
Put[expr,"temp.m"]; (* export the expression in a temporary file *)
Get["!"<>OptionValue[Path]<>"mjollnir --full-simplify temp.m" ](*execute Mjollnir*),
Resolve[expr,Reals] (* execute resolve *)
];

Test

By default, quantifier elimination is performed using Mjollnir. Run the following test to check if Mjollnir works correctly.

expr=ForAll[z,Implies[x≥0,And[Implies[x≥1,4≤x&&x≤9],Implies[x≤0,5≤z&&z≤10]]]&&Implies[x≤-1,And[Implies[x≥1,5≤z&&z≤10],Implies[x≤0,3≤x&&x≤5]]]];
QE[expr]

In case the test does not work you can set the Method as Automatic and use Resolve

QE[expr,Method->Automatic]

Weighted Kripke Structures Utilities

ApplyWeights[PWKS[model_],vals_]:=
Module[{lst,f},
f[w_,{i_}]:= w/;OddQ[i];
f[w_,{i_}]:= "<"<>ToString@ReplaceAll[ToExpression@w,vals]<>">"/;EvenQ[i];
lst= StringJoin@MapIndexed[f,StringSplit[model,{"<",">"}]]
];

The following function displays models. The function requires one to install Graphviz and the PVTool2. Once installed you shall set the installation directories as done below.

Options[ShowModel]={PVToolPath->"/Users/giovbacci/AAU/Uncertain-Quantitative-Systems/PVTool/",
GraphvizPath-> "/usr/local/Cellar/graphviz/2.40.1/bin/"};
ShowModel[PWKS[model_]]:=ShowModel[model];
ShowModel[WKS[model_]]:=ShowModel[model];
ShowModel[WUKS[model_,dcl_]]:=ShowModel[model];
ShowModel[model_String,OptionsPattern[]]:=
Module[{cmd = "!./pwks2wktool",out},
SetDirectory[OptionValue[PVToolPath]];
cmd = cmd <>" \""<>model<>"\" > model.dot";
ReadList[cmd,String];
Get["!"<>OptionValue[GraphvizPath]<>"dot -T pdf model.dot -o model.pdf"]; 
ResetDirectory[];
Show[Import[OptionValue[PVToolPath]<>"model.pdf","PDF"][[1]],ImageSize->Medium]
];

Sample WKSs from a given WUKS. The models will be saved in different files within the selected destination directory path

Options[SampleModels]={Error-> 0.005,Confidence->0.0001};
SampleModels[P:WUKS[model_,dcl_],path_,OptionsPattern[]]:=
Block[{data, n=SamplesNumber[OptionValue[Error],OptionValue[Confidence]],f,
X,Distr,adm = AdmissibleValuations[PWKS[model]]},
f[vals_,{index_}]:=
Module[{mapvals=Thread[Rule[X,vals]] },
If[adm/.mapvals,
Export["model"<>ToString[index]<>".wks",ApplyWeights[PWKS[model],mapvals],"String"],
Sequence[]
]];
CreateDirectory[path];
SetDirectory[path];
{X,Distr}=Transpose[dcl];
{Model -> P, Error->OptionValue[Error],Confidence->OptionValue[Confidence], Samples-> n,
Parameters-> X,Distributions-> Distr}>>"info.txt";
data=Round[10^6 *Transpose[RandomVariate[#,n]&/@Distr]];
MapIndexed[f,data];
ResetDirectory[];
"Done"
];

Model checking Weighted Kripke Structures under Uncertainties

Model Checking Parametric Weighted Kripke Structures

Interpret the output of the tool

Parse the output from the PVtool

ParseOutput[str_]:=ToExpression@StringReplace[str,{"{"-> "[","}"-> "]","INFTY"-> "Infinity"}];

Partial Evaluation of the parsed output

PVMAX[]:= 0;
PVMAX[a___,Infinity,b___]:=Infinity;
PVMAX[a___,0,b___]:= PVMAX[a,b]; 
PVMAX[x_]:=x; 

PVMIN[]:= Infinity;
PVMIN[a___,Infinity,b___]:=PVMIN[a,b];
PVMIN[a___,0,b___]:= 0;
PVMIN[x_]:=x;

Eval[NOT[a_ ≤ n_]]:=If[Not[Eval[a]≤Eval[n]],0,Infinity];
Eval[a_ ≤ n_]:= If[Eval[a]≤Eval[n],0,Infinity];
Eval[True]=0; Eval[False]:=Infinity;
Eval[x_PVMAX]:=Map[Eval,x,{1}];
Eval[x_PVMIN]:=Map[Eval,x,{1}];
Eval[x_]:=x;

Encode the numerical formula as an first-order logic formula in the existential (linear) real theory of the reals.

ToConstr[x_PVMAX≤q_]:=And@@Map[ToConstr[#≤q]&,x,{1}];
ToConstr[x_PVMAX>q_]:=And@@Map[ToConstr[#>q]&,x,{1}];
ToConstr[x_PVMIN≤q_]:=Or@@Map[ToConstr[#≤q]&,x,{1}];
ToConstr[x_PVMIN>q_]:=Or@@Map[ToConstr[#>q]&,x,{1}];
ToConstr[If[cond_,_,_]≤q_]:=ToConstr[cond];
ToConstr[If[cond_,_,_]>q_]:=ToConstr[Not[cond]];
ToConstr[x_Plus≤ q_]:=ToConstr[Distribute[x,PVMAX,Plus]≤q]/;Or@@Map[MatchQ[#,_PVMAX]&,x,{1}];
ToConstr[x_Plus≤ q_]:=ToConstr[Distribute[x,PVMIN,Plus]≤q]/;Or@@Map[MatchQ[#,_PVMIN]&,x,{1}];
ToConstr[x_]:=x;

Model Checking pWKS against bounded-WCTL

Compute the least fixed point of the EPDG associated with the pWKS P, state s, and the bounded WCTL formula phi. The function requires one to have installed Graphviz and the PVTool2. Once installed you shall set the installation directories as done below.

Options[PVtool]=Join[Options[ShowModel],{Evaluated->True, EPDG->False}];
PVtool[s_,P:PWKS[model_],phi_,OptionsPattern[]]:=
Module[{pvtoolcmd = "!./pvtool",out,dot},
SetDirectory[OptionValue[PVToolPath]];
pvtoolcmd = pvtoolcmd <>" \""<>model<>"\" \""<>phi<>"\" \""<>s<>"\"";
out=ReadList[pvtoolcmd,String]; (*execute pvtool*)
If[OptionValue[EPDG],(* show the Generated EPDG *)
Get["!"<>OptionValue[GraphvizPath]<>"dot -T pdf epdg.dot -o epdg.pdf"]; (* this does not always work !*)
Print@Show[Import[OptionValue[PVToolPath]<>"epdg.pdf","PDF"][[1]],ImageSize->Large];
];
ResetDirectory[];
If[OptionValue[Evaluated],
Eval[ParseOutput[Last[out]]],
Print@@out;
Last[out]
]
];

Symbolic representation of the set of evaluations solving the model checking problem

Options[Models]={Error-> 0.005,Confidence->0.0001,Verbose->False};
Models[P:PWKS[_],s_,phi_,OptionsPattern[]]:=
Simplify[Rationalize[ToConstr[PVtool[s,P,phi]≤0]&& AdmissibleValuations[P]]];

Symbolic representation of the admissible valuations

AdmissibleValuations[PWKS[model_]]:=
Module[{lst},
lst= StringSplit[model,{"<",">"}];
Simplify[And@@((ToExpression[#]≥0&)/@Extract[lst,Table[{2i},{i,Floor@Length[lst]/2}]])]
];

Example

P=PWKS["s2={a}<p> s2 s1={b}<3p+q>s2"]; (* the model *)
phi="Not E ((A X[10](a || b)) U[10] a)";

ShowModel[P]

PVtool["s1",P,phi,EPDG->True,Evaluated->True](* compute the least fixed point and the corresponding EPDG *)

Models[P,"s1",phi] (* compute the set valuations solving the model checking problem *)

Probabilistic Model Checking of Weight-Uncertain Kripke Structures

Monte Carlo Estimation

Number of samples based on the Chernoff bound

SamplesNumber[err_,conf_]:=Ceiling[-Log[conf]/(2 err^2)];

Estimate the event by sampling the random variables

Options[EstimateProb]={Error-> 0.005,Confidence->0.0001,Verbose->False};
EstimateProb[X_,Distr_, {c1_,c2_},OptionsPattern[]]:=
Block[{data, n=SamplesNumber[OptionValue[Error],OptionValue[Confidence]],pc1,pc2},
data=Transpose[RandomVariate[#,n]&/@Distr];
{pc1,pc2}=Mean/@Thread[Boole[{c1,c2}/.Thread[Rule[X,#]]]&/@data];
If[pc2>0,pc1/pc2,Indeterminate]
];
EstimateProb[X_,Distr_, constraint_,OptionsPattern[]]:=
Block[{data, n=SamplesNumber[OptionValue[Error],OptionValue[Confidence]]},
data=Transpose[RandomVariate[#,n]&/@Distr];
Mean[Boole[constraint/.Thread[Rule[X,#]]]&/@data]
];

Model Checking WUKS against bounded-WCTL

Models[J:WUKS[model_,dcl_],s_,phi_,opt:OptionsPattern[]]:=
Module[{csat,cadm=AdmissibleValuations[PWKS[model]],X,distr,tp,te,val,
err=OptionValue[Error],conf=OptionValue[Confidence]},
{tp,csat}=Timing[Models[PWKS[model],s,phi]];
Switch[csat,
True,{te,val}={0,1},
False,{te,val}={0,0}, 
_,
{X,distr}=Thread[dcl];
{te,val}=Timing[EstimateProb[X,distr,{csat,cadm},opt]];
];
If[OptionValue[Verbose],
Print["Event:  E = (",csat,") [computed in ", tp,"sec]"];
Print["Admissible Valuations: V = (",cadm,")"];
Print["Samples: ", SamplesNumber[err,conf], ", Error: ", err,", Confidence: ", 100(1-conf), "%"];
Print["Estimation of P[X∈E|X∈V] [computed in ", te,"sec]"]
];
val
];

Example

dcl = {{p,UniformDistribution[{0,1}]},{q,UniformDistribution[{0,1}]}};
J=WUKS["s2={a} s1={b}<3p+q>s2",dcl];
phi="A ((E X[10](a || b)) U[2] a)";

Models[J,"s1",phi]

Examples

Lawn Mower Example

Consider the PWKS below representing a grass field with different routes a lawn mower can take from the starting state s0 to s6 where the grass can be dumped. The weights on the transitions represent the amount of grass that is accumulated in the container when selecting a particular route.

lawnMower=PWKS["s0={mow}<p>s1 <p>s2 <p>s3 
s1={mow} <q>s4 
s2={mow} <2q>s4 
s3={mow} <q>s5 
s4={mow} <r>s5 <2q>s6 
s5={mow} <p>s6 
s6={dump} <0>s6"];

ShowModel[lawnMower]

Assume that the lawn mower breaks when it is forced to store more than 6.5 units of grass, then the property “the grass is always dumped before the lawn mower breaks, irrelevant of the selected route” is expressed by the WCTL formula ϕ = ∀(mow (Subscript[U, ≤6.5]) dump).

ϕ="A (mow U[6.5] dump)";

To model check the PWKS lawnMower against ϕ, we can invoke the PVTool. The PVTool will construct an Extended Parametric Dependency Graph (EPGD) and symbolically compute its least fixed point assignment associated to the root configuration ⟨s0, ϕ⟩. 
Optionally, one can visualise the EPDG by setting the “EPDG” option to True (by default the option is disabled).

PVtool["s0",lawnMower,ϕ,EPDG->True]

The set of admissible valuations for the lawnMower model can be retrieved as follows

AdmissibleValuations[lawnMower]

The set of valuations satisfying ϕ, i.e., [[lawnMower, s0 ⊨ ϕ]], can be represented by the following linear constraints

Models[lawnMower,"s0",ϕ]

Assume that we have measured p ≃ 2 ± ϵ, q ≃ 1 ± ϵ and r ≃ 0 ± ϵ where ϵ  represents the the absolute measurement error of our parameters. One can determine if the model is robust with respect to ϕ by checking that all possible measurements values lay in [[lawnMower, s0 ⊨ ϕ]]. 
This can be expressed by the following first-order formula in the theory of real linear arithmetic

robust =ForAll[{p,q,r},
Implies[0<=p &&p ≤2+err && 0≤q && q ≤1+err && 0<=r &&r≤err,
Models[lawnMower,"s0",ϕ]]]

we can eliminate the quantifiers and obtain an equivalent formula

QE[robust]

The above formula indicates that to ensure robustness, the error ϵ  shall not exceed 1/10.

Assume now that the measurements for p, q and, r are associated to a distribution, i.e., they are modelled as real-valued random variables. The lawnMower shall be now modelled as a Weight-Uncertain Kripke Structure (WUKS). 

Consider for instance that p~N(2,ϵ), q~Unif(1-ϵ,1+ϵ) and, r~N(0,ϵ) where ϵ= 1/10. We can declare the corresponding WUKS, called uncLawnMower, as follows

dcl={{p,NormalDistribution[2,0.1]},{q,UniformDistribution[{0.9,1.1}]},{r,NormalDistribution[0,0.1]}};
uncLawnMower=WUKS[lawnMower[[1]],dcl];

Note that the tool supports any of the distributions provided by Mathematica (cf. Random Variables).

Now we can estimate the probability that uncLawnMower satisfies ϕ as follows. In this example we have set the error and the confidence for the estimation respectively to 0.005 and 0.001 by setting the options Error and Confidence (default values are Error -> 0.005 Confidence -> 0.0001). Futhermore one can set the option Verbose to True to printout some information regarding the estimation.

N@Models[uncLawnMower,"s0",ϕ,Verbose->True,Error->0.005, Confidence->0.001]

Factory Example

factorymodel ="
factory = {}<10>factory <2>s1 <10>s2 <5>s2
s1 ={} <10>s1 <10>s3 <5>s3
s2 ={} <10>s2 <2>s3 <2>s4 <2>s5
s3 ={} <10>s3 <2>s5
s4 ={sad,happy} <10>s4 <2>s5
s5 ={sad,happy} <10>s5 <0>s6 <0>s9
s6 ={goodproduct} <10>s6 <10>s7 <5>s7
s7 ={goodproduct} <10>s7 <2>s8
s8 ={goodproduct,sad,happy} <10>s8
s9 ={badproduct} <10>s9 <10>s10 <5>s10
s10 ={badproduct} <10>s10 <2>s11
s11 ={badproduct,sad,happy} <10>s11";

psi1="E (true U[18] happy)";
psi2="E (true U[50] badproduct)";

PVtool["factory",PWKS[factorymodel],psi2,EPDG->True,Evaluated->False]

Models[PWKS[factorymodel],"factory",psi1]