Using multivalue DDs to solve multistate reliability quantification

@DCTLib, do you recall this discussion below? You suggested a recursive equation, which was the right approach.

Cudd_PrintMinterm, accessing the individual minterms in the sum of products

Now, I am considering multistate reliability, where we can have either not fail or fail to n-1 different states, with n >= 2. Tulip-dd implements MDDs as described in:

https://github.com/tulip-control/dd/blob/master/doc.md#multi-valued-decision-diagrams-mdd

https://github.com/tulip-control/dd/issues/71

https://github.com/tulip-control/dd/issues/66

In the diagrams in the drawings below, we have defined an MDD declared by:

aut.declare_variable(x=(0,3)) u = aut.add_expr(‘x=i’)

Each value/state of the multi-value variable (MSV) x, x=0, x=1, x=2, or x=3 leads to a specific BDD as shown in the diagrams at the bottom, taking a four-state variable x as example here. The notation is that state 0 represents the normal state and x can fail to different states 1, 2, and 3. The failure probabilities are assigned in table below. In the BDDs below, we (and tulip as well) use the binary coding with two bits x_1 and x_0 to represent each state/value of the MSV. The least significant bit (LSB), i.e., x_0, is always the ancestor. Each of the BDD diagrams below is a representation of a specific value, or state.

To quantify the BDD of a specific state, i.e., the top node, we must know probabilities of binary variables x_0 and x_1 taking different branches (then or else) in the BDD. These branch probabilities are not given directly but need to be calculated according to the BDD structure.

The key here is that the child node probabilities and the branch probabilities of the parent node must be known prior to the calculation of the parent node probability. In the previous BDD quantification, we knew the probabilities of branches from node x_1 to leaf nodes when calculating node x_1 probability. We did not need to know how node x_1 was connected to node x_0. Now, for this four-state variable x, we need to know how node x_1 is connected to node x_0, the binary variable representing the least significant bit, to determine the probabilities of branches from node x_1 to leaf nodes. The question is how to implement this?

Here’s one attempt that falls short:

import numpy as np
from omega.symbolic import temporal as trl
aut = trl.Automaton()
# Example of function that returns branch probabilities
def prntr(v, pvars):
    
    assert sum(pvars)==1.0,"Probs must add to 1!"

    if (v.high).var == 'x_1':
        if (v.low) == aut.true:
            return pvars[1] + pvars[3], pvars[1]
        else:
            return pvars[1] + pvars[3], pvars[3]
    if (v.low).var == 'x_1':
        if (v.low).negated:
            return pvars[1] + pvars[3], pvars[0]
        else:
            return pvars[1] + pvars[3], pvars[2]
aut.declare_variables(x=(0, 3))
u=aut.add_expr('x = 3')
pvars = np.array([0.1, 0.2, 0.3, 0.4])
prntr(u,pvars)

The key here is that the child node probabilities and the branch probabilities of the parent node must be known prior to the calculation of the parent node probability.

Yes, exactly. In this case, a fully recursive bottom-up computation, like normally done with BDDs, will not work for the reason that you wrote.

However, the approach will start to work again when you treat the variables that together form a state to be a block. So in your recursive function for the probability calculation, whenever you encounter a variable for a block, you treat the node and the successor nodes for the same state component as a block and only recurse when you encounter a node not belonging to the block.

Note that this approach requires that the variables for the state appear continuously in the variable ordering. For the CUDD library, you can constrain the automatic variable reordering to guarantee this.

The following code is a modification of yours implementing this idea:

#!/usr/bin/env python3
import numpy as np
from omega.symbolic import temporal as trl
aut = trl.Automaton()

# Example of function that returns branch probabilities
# Does not yet use result caching and does not yet support assigning probabilities 
# to more than one state variable set
def prntr(v, pvars):
    
    assert abs(sum(pvars)-1.0)<=0.0001,"Probs must add to 1!"
    
    if v == aut.true:
        return 1.0
    elif v == aut.false:
        return 0.0
        
    if v.var in ["x_0","x_1"]:
        thisSum = 0.0
        # Compute probabilities
        for i,p in enumerate(pvars):
            
            # Find successor
            # Assuming that x_2 always comes after x_1
            thisV = v
            negate = thisV.negated
            if thisV.var == 'x_0':
                if i & 1:
                    thisV = thisV.high
                else:
                    thisV = thisV.low
                    negate = negate ^ thisV.negated
            if thisV.var == 'x_1':
                if i & 2:
                    thisV = thisV.high
                else:
                    thisV = thisV.low
            if negate:
                thisSum += p*prntr(~thisV, pvars)
            else:
                thisSum += p*prntr(thisV, pvars)
        
        return thisSum

    # TODO: What is the semantics for variables outside of the current block?
    return 0.5*prntr(v.high, pvars) + 0.5*prntr(v.low, pvars)

pvars = np.array([0.1, 0.2, 0.3, 0.4])
aut.declare_variables(x=(0, 3))
u= aut.add_expr('x = 0')
print(prntr(u,pvars))
u2 = aut.add_expr('x = 3') | aut.add_expr('x = 2')
print(prntr(u2,pvars))