Next Steps

Up until now, we showcased various analysis approaches for different model types in the previous steps. This notebook presents typical use cases for the Python API: Using symbolic model representations, using preprocessing algorithms, configuring solvers, and implementing own prototypes on top of existing routines.

Preparation

As usual, we load the Prism Orchard model.

[1]:

# Import prism models from the previous step
import stormvogel
import stormpy
from orchard.orchard_builder import build_prism

[2]:

orchard_prism, prism_program = build_prism()

Model representation

So far, we have represented the MDP through an explicit (sparse) transition matrix. Instead, we can build the MDP symbolically and represent it through decision diagrams (DD).

[3]:

orchard_symbolic = stormpy.build_symbolic_model(prism_program)
print(orchard_symbolic)

--------------------------------------------------------------
Model type:     MDP (symbolic)
States:         22469 (85 nodes)
Transitions:    44954 (916 nodes)
Choices:        29354
Reward Models:  rounds
Variables:      rows: 7 meta variables (19 DD variables), columns: 7 meta variables (19 DD variables), nondeterminism: 12 meta variables (12 DD variables)
Labels:         6
   * deadlock -> 0 state(s) (1 nodes)
   * init -> 1 state(s) (20 nodes)
   * AllCherriesPicked
   * PlayersWon
   * RavenOneAway
   * RavenWon
--------------------------------------------------------------

The symbolic MDP has the same size as the sparse MDP, but instead of explicitly storing nearly 45,000 transitions, the function is represented through a DD with 946 nodes.

We can perform model checking on this symbolic representation and obtain the same winning probability.

[4]:

formula = stormpy.parse_properties('Pmax=? [F "PlayersWon"]')[0]
symbolic_result = stormpy.model_checking(
    orchard_symbolic, formula, only_initial_states=True
)
filter = stormpy.create_filter_initial_states_symbolic(orchard_symbolic)
symbolic_result.filter(filter)
print("Maximal probability: {}".format(symbolic_result))

WARN  (SymbolicMinMaxLinearEquationSolver.cpp:64): Selected method is not supported for this solver, switching to value iteration.
WARN  (SymbolicMinMaxLinearEquationSolver.cpp:64): Selected method is not supported for this solver, switching to value iteration.
Maximal probability: 0.631357

Bisimulation

Bisimulation minimization can be applied to minimize the state space without affecting the model behaviour.

[5]:

formula = stormpy.parse_properties('Pmax=? [F "PlayersWon"]')
print(
    "Model with {} states and {} transitions".format(
        orchard_prism.nr_states, orchard_prism.nr_transitions
    )
)
orchard_bisim = stormpy.perform_bisimulation(
    orchard_prism, formula, stormpy.BisimulationType.STRONG
)
print(
    "Model with {} states and {} transitions".format(
        orchard_bisim.nr_states, orchard_bisim.nr_transitions
    )
)

Model with 22469 states and 44954 transitions
Model with 956 states and 2446 transitions

For the orchard model, we can reduce the state space size by \(95\%\) from \(22,469\) states to \(956\) states. This reduction also speeds up subsequent analysis on the reduced model.

Model checking algorithms

The underlying model checking algorithms in Storm are configured through so-called environments. For example, we can set the precision requirement from the default value of \(10^{−6}\) to \(0.1\) and see how that affects the resulting probability.

[6]:

# Change precision
env = stormpy.Environment()
prec = stormpy.Rational(0.1)
env.solver_environment.minmax_solver_environment.precision = prec
result = stormpy.model_checking(orchard_prism, formula[0], environment=env)
print(result.at(orchard_prism.initial_states[0]))

0.5815061686029692

We can see that the resulting probability \(0.5815\) is within the precision guarantee of \(0.1\).

We can also change the underlying algorithm and compare different approaches such as value iteration, policy iteration and optimistic value iteration.

[7]:

# Change algorithm
import time

methods = [
    stormpy.MinMaxMethod.value_iteration,
    stormpy.MinMaxMethod.policy_iteration,
    stormpy.MinMaxMethod.optimistic_value_iteration,
]
for m in methods:
    env = stormpy.Environment()
    env.solver_environment.minmax_solver_environment.method = m
    start = time.time()
    result = stormpy.model_checking(
        orchard_prism, formula[0], environment=env, extract_scheduler=True
    )
    print(f"Method: {m}")
    print(f"Result: {result.at(orchard_prism.initial_states[0])}")
    print(f"Time: {time.time() - start:.2}s")

Method: MinMaxMethod.value_iteration
Result: 0.6313572986959963
Time: 0.066s
Method: MinMaxMethod.policy_iteration
Result: 0.6313573066006365
Time: 0.3s
Method: MinMaxMethod.optimistic_value_iteration
Result: 0.6313576699315778
Time: 0.036s

We see that all methods provide the same result, but their timings differ. For example, policy iteration is typically slower.

Running the above cell multiple times allows to get better insights into the performance of the different algorithms.

Writing an LP-based MDP Model Checker

One further step is extending stormpy with own algorithms. Extending stormpy allows to profit from existing data structure and efficient algorithms when integrating own ideas. We exemplify this with a fundamental approach to solve MDPs using a linear program (LP) solver, an approach that is popular due to the flexibility to add additional constraints.

The LP encoding below computes reachability probabilities.

Given a target label \(B \in AP\), this LP characterizes the variables \((x_s)_{s \in S}\) as follows:

Minimize \(\sum_{s \in S} x_s\) such that:

  1. If \(B \in L(s)\), then \(x_s =1\),

  2. If there is no path from \(s\) to any state with label \(B\), then \(x_s=0\),

  3. Otherwise, \(0 \leq x_s \leq 1\), and for all actions \(\alpha \in \textit{Act}(s)\): \(x_s \geq \sum_{s' \in S} \mathbb{P}(s, \alpha, s') \cdot x_{s'}\)

We encode the LP by using functionality of stormpy. We first compute all states that satisfy the requirements in (1) and (2) using graph algorithms provided in stormpy. For instance, we compute the states \(s\) which have to path to label \(B\) using the function compute_prob01max_states(). We then construct and solve the LP by inspecting the model’s transition matrix.

[8]:

# These are the target states
players_won = stormpy.parse_properties('"PlayersWon"')[0].raw_formula
psi_states = stormpy.model_checking(orchard_prism, players_won).get_truth_values()

# These are the states that can never reach the target states
phi_states = stormpy.BitVector(orchard_prism.nr_states, True)
prob0max_states = stormpy.compute_prob01max_states(
    orchard_prism, phi_states, psi_states
)[0]

# SCIP is an LP solver
from pyscipopt import Model

m = Model()

# Create a variable for each state
num_states = orchard_prism.nr_states
state_vars = [m.addVar(f"x_{i}", lb=0, ub=1) for i in range(num_states)]

# Encode LP
for state in range(num_states):
    if psi_states.get(state):  # Case 1
        m.addCons(state_vars[state] == 1)
    elif prob0max_states.get(state):  # Case 2
        m.addCons(state_vars[state] == 0)
    else:  # Case 3
        for row in orchard_prism.transition_matrix.get_rows_for_group(state):
            summed_prob = 0
            for transition in orchard_prism.transition_matrix.get_row(row):
                prob = transition.value()
                next_state = transition.column
                summed_prob += prob * state_vars[next_state]
            m.addCons(state_vars[state] >= summed_prob)

Lastly, we can solve the LP and obtain the same winning probability as before. We also compare the LP result with the model checking result.

[9]:

# Solve LP
m.setObjective(sum(state_vars), sense="minimize")
m.optimize()
sol = m.getBestSol()
result_lp = sol[state_vars[orchard_prism.initial_states[0]]]

# Compare with default VI
properties = stormpy.parse_properties_without_context('Pmax=? [F "PlayersWon"]')
vi_result = stormpy.model_checking(orchard_prism, properties[0].raw_formula)
result_vi = vi_result.at(orchard_prism.initial_states[0])

print(result_lp, result_vi)
assert abs(result_lp - result_vi) < 1e-6  # => True

presolving:
(round 1, fast)       16318 del vars, 16318 del conss, 0 add conss, 1253 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 2, fast)       16347 del vars, 16351 del conss, 0 add conss, 1471 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 3, fast)       16383 del vars, 16487 del conss, 0 add conss, 2510 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 4, fast)       16496 del vars, 16679 del conss, 0 add conss, 5251 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 5, fast)       16624 del vars, 16979 del conss, 0 add conss, 9886 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 6, fast)       16868 del vars, 17410 del conss, 0 add conss, 14995 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 7, fast)       17133 del vars, 18067 del conss, 0 add conss, 20025 c0.6313573066006342 0.6313566764945008
hg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 8, fast)       17550 del vars, 18959 del conss, 0 add conss, 24516 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 9, fast)       18047 del vars, 20139 del conss, 0 add conss, 28476 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 10, fast)       18674 del vars, 21479 del conss, 0 add conss, 31530 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 11, fast)       19377 del vars, 23106 del conss, 0 add conss, 34095 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 12, fast)       20099 del vars, 24631 del conss, 0 add conss, 36007 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints
(round 13, fast)       20784 del vars, 26084 del conss, 0 add conss, 37365 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clq