# Fault Tree with qubovert

*qubovert* must be pip installed.

We will encode the Fault Tree problem shown below. The inputs are `e0`, `e1`, `e2`, and `e3`, which are boolean variables; for example, `e0 == 1` if Error 0 occurs, otherwise it is 0. The output is `p0`, which is a boolean indicating top level failure; it is 1 if a top level failure occurs, otherwise it is 0. The goal of this analysis is to find the minimum number of errors (`e0`, `e1`, `e2`, `e3`) that must be 1 in order for `p0` to be 1.



The Fault Tree image above can be summarized by the following:

\begin{align*}
z0 &= e0 \text{ OR }  e3\\
z1 &= e2 \text{ OR } z0\\
z2 &= e0 \text{ OR } e1\\
p0 &= z1 \text{ AND } z2
\end{align*}

We will solve this problem in two ways:

1. <a href="#sat">Solving with the `qubovert.sat` library</a>
1. <a href="#pcbo">Solving with `qubovert.PCBO`</a>

<div id="sat" />

## Solving with the `qubovert.sat` library

Let's create the variables and use the `qubovert.sat` libary to encode `p0`.


```python
from qubovert import boolean_var
from qubovert.sat import OR, AND

e0, e1 = boolean_var('e0'), boolean_var('e1')
e2, e3 = boolean_var('e2'), boolean_var('e3')

z0 = OR(e0, e3)
z1 = OR(z0, e2)
z2 = OR(e0, e1)
p0 = AND(z1, z2)

print(p0)
```

    {('e0',): 1, ('e0', 'e1', 'e3'): -1, ('e1', 'e3'): 1, ('e0', 'e1', 'e3', 'e2'): 1, ('e1', 'e3', 'e2'): -1, ('e0', 'e1', 'e2'): -1, ('e1', 'e2'): 1}


We want to find the minimum number of errors that will lead to a top level error. Thus, we want to minimize `H`.


```python
H = e0 + e1 + e2 + e3

print(H)
```

    {('e0',): 1, ('e1',): 1, ('e2',): 1, ('e3',): 1}


We now subject `H` to the constraint that `p0 == 1`, or equivalently `1 - p0 == 0`, where we notice that `1 - p0` is bounded below by 0 and above by 1. We will add a penalty `lam` to the PCBO to enforce this constraint. For now, let's make it a symbol that we can tune later.


```python
#!pip install sympy
from sympy import Symbol

lam = Symbol("lam", positive=True)

H.add_constraint_eq_zero(1 - p0, lam=lam, bounds=(0, 1))

print("PCBO:\n", H, "\n")
print("Constraints:\n", H.constraints)
```

    PCBO:
     {('e0',): 1 - lam, ('e1',): 1, ('e2',): 1, ('e3',): 1, ('e0', 'e1', 'e3'): lam, ('e1', 'e3'): -lam, ('e0', 'e1', 'e3', 'e2'): -lam, ('e1', 'e3', 'e2'): lam, ('e0', 'e1', 'e2'): lam, ('e1', 'e2'): -lam, (): lam} 
    
    Constraints:
     {'eq': [{('e0',): -1, ('e0', 'e1', 'e3'): 1, ('e1', 'e3'): -1, ('e0', 'e1', 'e3', 'e2'): -1, ('e1', 'e3', 'e2'): 1, ('e0', 'e1', 'e2'): 1, ('e1', 'e2'): -1, (): 1}]}


Notice there is one equality constraint coming from the requirement that `1 == p0`. The `H.is_solution_valid` function will take in a proposed solution to the problem and ensure that this constraint is satisfied.

The `'eq'` key of the constraints dictionary indicates that the quantity equals zero, and the `'lt'` key of the constraints dictionary indicates that the quantity is less than zero. Other possible keys are `'le'`, `'gt'`, and `'ge'`. See the docstrings for `PCBO.add_constraint_eq_zero`, `PCBO.add_constraint_lt_zero`, `PCBO.add_constraint_le_zero`, `PCBO.add_constraint_gt_zero`, and `PCBO.add_constraint_ge_zero` for info.

For testing purposes, let's solve this bruteforce to make sure everything is working.


```python
solutions = H.solve_bruteforce(all_solutions=True)
print(solutions)
```

    [{'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0}]


Notice that there is one unique solution that minimizes the objective function and obeys all the constraints. If just the error `e0` occurs, then a top level failure will occur.


```python
solution = solutions[0]
print("Minimum number of failures that leads to a top level failure:", H.value(solution))
print("p0 =", p0.value(solution))
```

    Minimum number of failures that leads to a top level failure: 1
    p0 = 1


Now let's solve this problem with a generic QUBO solver. Notice that the degree of problem is more than two, making `H` not a natural Quadratic Unconstrained Boolean Optimization Problem (QUBO).


```python
H.degree
```




    4



We can convert it to a QUBO (note that there are some options for the reduction from PUBO to QUBO, see the `H.to_qubo` method for details). Ancilla bits will need to be added, and bit labels are mapped to integers.


```python
Q = H.to_qubo()
print("num PUBO variables", H.num_binary_variables)
print("num QUBO variables", Q.num_binary_variables)
print()
print(Q)
```

    num PUBO variables 4
    num QUBO variables 6
    
    {(0,): 1 - lam, (1,): 1, (2,): 1, (3,): 1, (4,): 9*lam + 9, (1, 3): 2*lam + 3, (1, 4): -6*lam - 6, (3, 4): -6*lam - 6, (0, 4): lam, (5,): 6*lam + 6, (0, 2): 2*lam + 2, (0, 5): -4*lam - 4, (2, 5): -4*lam - 4, (4, 5): -lam, (2, 4): lam, (1, 5): lam, (1, 2): -lam, (): lam}


For testing purposes, let's solve this with bruteforce to see what the proper value of $\lambda$ should be to enforce the constraints. Notice how we remap the QUBO solution to the PCBO solution with `H.convert_solution(x)`.


```python
for l in (1, 2, 3):
    Q_temp = Q.subs({lam: l})
    solutions = Q_temp.solve_bruteforce(all_solutions=True)
    solutions = [H.convert_solution(x) for x in solutions]
    print('lam', l)
    for s in solutions:
        print("\t", s, "is", "valid" if H.is_solution_valid(s) else "invalid")
    print()
```

    lam 1
    	 {'e0': 0, 'e1': 0, 'e2': 0, 'e3': 0} is invalid
    	 {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0} is valid
    
    lam 2
    	 {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0} is valid
    
    lam 3
    	 {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0} is valid
    


We see that $\lambda = 2$ is sufficient to enforce the constraints. So let's update our QUBO.


```python
Q_good = Q.subs({lam: 2})
```

Now let's solve the QUBO with D'Wave's simulated annealer.


```python
#!pip install dwave-neal
from neal import SimulatedAnnealingSampler

sampler = SimulatedAnnealingSampler()
```

Note that their software package takes in a specific form for QUBOs, namely, the keys of the dictionary must be two element tuples. This form can be accessed from `Q` and `Q_good` with `Q.Q` or `Q_good.Q`.


```python
qubo_sample = sampler.sample_qubo(Q_good.Q, num_reads=100)
print("objective function:", qubo_sample.first.energy + Q_good.offset, "\n")

qubo_solution = qubo_sample.first.sample
print("qubo solution:", qubo_solution, "\n")

solution = H.convert_solution(qubo_solution)
print("pcbo solution:", solution)
print("objective function:", H.value(solution), "\n")

print("The solution is", "valid" if H.is_solution_valid(solution) else "invalid")
```

    objective function: 1.0 
    
    qubo solution: {0: 1, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0} 
    
    pcbo solution: {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0}
    objective function: 1 
    
    The solution is valid


This matches the result of `H.solve_bruteforce()`. Recall that the objective function is equal to the minimum number of failures that will lead to a top level failure.

Now we'll solve an QUSO formulation of our problem. Again we'll take $\lambda = 2$.


```python
L = H.to_quso().subs({lam: 2})
# note that we cannot do H.subs({lam: 2}).to_quso()!! This is because H.subs({lam: 2})
# creates a new PCBO object, and it's mapping from variables labels to integers may be
# different than H's mapping. For example, try H.mapping == H.subs({lam: 2}).mapping a
# few times. They will often be different.
print("num PUBO variables", H.num_binary_variables)
print("num QUSO variables", L.num_binary_variables)
print()
print(L)
```

    num PUBO variables 4
    num QUSO variables 6
    
    {(0,): 1.5, (): 14.25, (1,): 2.25, (2,): 1.0, (3,): 2.25, (4,): -5.0, (1, 3): 1.75, (1, 4): -4.5, (3, 4): -4.5, (0, 4): 0.5, (5,): -3.0, (0, 2): 1.5, (0, 5): -3.0, (2, 5): -3.0, (4, 5): -0.5, (2, 4): 0.5, (1, 5): 0.5, (1, 2): -0.5}


Similar to their QUBO solver, D'Wave's QUSO solver accepts a specific form for QUSO models, namely a linear term dictionary and a quadratic term dictionary. These can be accessed with `L.h` and `L.J`.


```python
quso_sample = sampler.sample_ising(L.h, L.J, num_reads=100)
print("objective function:", quso_sample.first.energy + L.offset, "\n")

quso_solution = quso_sample.first.sample
print("quso solution:", quso_solution, "\n")

solution = H.convert_solution(quso_solution)
print("pcbo solution:", solution)
print("objective function:", H.value(solution), "\n")

print("The solution is", "valid" if H.is_solution_valid(solution) else "invalid")
```

    objective function: 1.0 
    
    quso solution: {0: -1, 1: 1, 2: 1, 3: 1, 4: 1, 5: 1} 
    
    pcbo solution: {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0}
    objective function: 1 
    
    The solution is valid


Again this matches the result of `H.solve_bruteforce()`.

<div id="pcbo" />

## Solving with `qubovert.PCBO`

We want to find the minimum number of errors that will lead to a top level error. Thus, we want to minimize `H`. We could make `H` the same way that we did in the previous section, by creating variables with ``qubovert.boolean_var``, but instead for illistration we will make it as a dictionary instead.


```python
from qubovert import PCBO

H = PCBO(
    {(x,): 1 for x in ('e0', 'e1', 'e2', 'e3')}
)

print(H)
```

    {('e0',): 1, ('e1',): 1, ('e2',): 1, ('e3',): 1}


Now let's enforce the constraints. We will need to enforce the constraints with a penalty factor $\lambda$. Let's create a symbol here that we can tune later.


```python
#!pip install sympy
from sympy import Symbol

lam = Symbol("lam", positive=True)
```

Now let's enforce the constraints (reproduced here for reference):
\begin{align*}
z0 &= e0 \text{ OR }  e3\\
z1 &= e2 \text{ OR } z0\\
z2 &= e0 \text{ OR } e1
\end{align*}


```python
H.add_constraint_eq_OR(
    'z0', 'e0', 'e3', lam=lam
).add_constraint_eq_OR(
    'z1', 'e2', 'z0', lam=lam
).add_constraint_eq_OR(
    'z2', 'e0', 'e1', lam=lam
)

print(H)
```

    {('e0',): 2*lam + 1, ('e1',): lam + 1, ('e2',): lam + 1, ('e3',): lam + 1, ('z0',): 2*lam, ('e0', 'e3'): lam, ('e0', 'z0'): -2*lam, ('z0', 'e3'): -2*lam, ('z1',): lam, ('z0', 'e2'): lam, ('e2', 'z1'): -2*lam, ('z0', 'z1'): -2*lam, ('z2',): lam, ('e0', 'e1'): lam, ('e0', 'z2'): -2*lam, ('e1', 'z2'): -2*lam}


Finally, we want to make $p0 = z1 \text{ AND } z2$ energetically favorable. We can do this with the following.


```python
H.add_constraint_AND('z1', 'z2', lam=lam)
print("PCBO:\n", H, "\n")
print("Constraints:\n", H.constraints)
```

    PCBO:
     {('e0',): 2*lam + 1, ('e1',): lam + 1, ('e2',): lam + 1, ('e3',): lam + 1, ('z0',): 2*lam, ('e0', 'e3'): lam, ('e0', 'z0'): -2*lam, ('z0', 'e3'): -2*lam, ('z1',): lam, ('z0', 'e2'): lam, ('e2', 'z1'): -2*lam, ('z0', 'z1'): -2*lam, ('z2',): lam, ('e0', 'e1'): lam, ('e0', 'z2'): -2*lam, ('e1', 'z2'): -2*lam, ('z2', 'z1'): -lam, (): lam} 
    
    Constraints:
     {'eq': [{('z0',): 1, ('e0',): 1, ('e3',): 1, ('e0', 'e3'): 1, ('e0', 'z0'): -2, ('z0', 'e3'): -2}, {('z1',): 1, ('e2',): 1, ('z0',): 1, ('z0', 'e2'): 1, ('e2', 'z1'): -2, ('z0', 'z1'): -2}, {('z2',): 1, ('e0',): 1, ('e1',): 1, ('e0', 'e1'): 1, ('e0', 'z2'): -2, ('e1', 'z2'): -2}, {('z2', 'z1'): -1, (): 1}]}


The `H.is_solution_valid` function will take in a proposed solution to the problem and ensure that these constraints are satisfied.

For testing purposes, let's solve this bruteforce to make sure everything is working.


```python
solutions = H.solve_bruteforce(all_solutions=True)
print(solutions)
```

    [{'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 1, 'z1': 1, 'z2': 1}]


Notice that there is one unique solution that minimizes the objective function and obeys all the constraints. If just the error `e0` occurs, then a top level failure will occur.


```python
solution = solutions[0]
print("Minimum number of failures that leads to a top level failure:", H.value(solution))
print("p0 =", H.value(solution))
```

    Minimum number of failures that leads to a top level failure: 1
    p0 = 1


Now let's solve this problem with a generic QUBO solver. Notice that the degree of problem is two, making `H` a natural Quadratic Unconstrained Boolean Optimization Problem (QUBO).


```python
H.degree
```




    2




```python
Q = H.to_qubo()
print(Q)
```

    {(0,): 2*lam + 1, (1,): lam + 1, (2,): lam + 1, (3,): lam + 1, (4,): 2*lam, (0, 3): lam, (0, 4): -2*lam, (3, 4): -2*lam, (5,): lam, (2, 4): lam, (2, 5): -2*lam, (4, 5): -2*lam, (6,): lam, (0, 1): lam, (0, 6): -2*lam, (1, 6): -2*lam, (5, 6): -lam, (): lam}


For testing purposes, let's solve this with bruteforce to see what the proper value of $\lambda$ should be to enforce the constraints. Notice how we remap the QUBO solution to the PCBO solution with `H.convert_solution(x)`.


```python
for l in (1, 2, 3):
    Q_temp = Q.subs({lam: l})
    solutions = Q_temp.solve_bruteforce(all_solutions=True)
    solutions = [H.convert_solution(x) for x in solutions]
    print('lam', l)
    for s in solutions:
        print("\t", s, "is", "valid" if H.is_solution_valid(s) else "invalid")
    print()
```

    lam 1
    	 {'e0': 0, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 0, 'z1': 0, 'z2': 0} is invalid
    	 {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 1, 'z1': 1, 'z2': 1} is valid
    
    lam 2
    	 {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 1, 'z1': 1, 'z2': 1} is valid
    
    lam 3
    	 {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 1, 'z1': 1, 'z2': 1} is valid
    


We see that $\lambda = 2$ is sufficient to enforce the constraints. So let's update our QUBO.


```python
Q_good = Q.subs({lam: 2})
```

Now let's solve the QUBO with D'Wave's simulated annealer.


```python
#!pip install dwave-neal
from neal import SimulatedAnnealingSampler

sampler = SimulatedAnnealingSampler()
```

Note that their software package takes in a specific form for QUBOs, namely, the keys of the dictionary must be two element tuples. This form can be accessed from `Q` and `Q_good` with `Q.Q` or `Q_good.Q`.


```python
qubo_sample = sampler.sample_qubo(Q_good.Q, num_reads=100)
print("objective function:", qubo_sample.first.energy + Q_good.offset, "\n")

qubo_solution = qubo_sample.first.sample
print("qubo solution:", qubo_solution, "\n")

solution = H.convert_solution(qubo_solution)
print("pcbo solution:", solution)
print("objective function:", H.value(solution), "\n")

print("The solution is", "valid" if H.is_solution_valid(solution) else "invalid")
```

    objective function: 1.0 
    
    qubo solution: {0: 1, 1: 0, 2: 0, 3: 0, 4: 1, 5: 1, 6: 1} 
    
    pcbo solution: {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 1, 'z1': 1, 'z2': 1}
    objective function: 1 
    
    The solution is valid


This matches the result of `H.solve_bruteforce()`. Recall that the objective function is equal to the minimum number of failures that will lead to a top level failure.

Now we'll solve an QUSO formulation of our problem. Again we'll take $\lambda = 2$.


```python
L = H.to_quso().subs({lam: 2})
# note that we cannot do H.subs({lam: 2}).to_quso()!! This is because H.subs({lam: 2})
# creates a new PCBO object, and it's mapping from variables labels to integers may be
# different than H's mapping. For example, try H.mapping == H.subs({lam: 2}).mapping a
# few times. They will often be different.
print("num QUSO variables", L.num_binary_variables)
print()
print(L)
```

    num QUSO variables 7
    
    {(0,): -1.5, (): 8.0, (1,): -1.0, (2,): -1.0, (3,): -1.0, (0, 3): 0.5, (0, 4): -1.0, (3, 4): -1.0, (2, 4): 0.5, (4,): 0.5, (2, 5): -1.0, (4, 5): -1.0, (5,): 1.5, (0, 1): 0.5, (0, 6): -1.0, (1, 6): -1.0, (6,): 1.5, (5, 6): -0.5}


Similar to their QUBO solver, D'Wave's QUSO solver accepts a specific form for QUSO models, namely a linear term dictionary and a quadratic term dictionary. These can be accessed with `L.h` and `L.J`.


```python
quso_sample = sampler.sample_ising(L.h, L.J, num_reads=100)
print("objective function:", quso_sample.first.energy + L.offset, "\n")

quso_solution = quso_sample.first.sample
print("quso solution:", quso_solution, "\n")

solution = H.convert_solution(quso_solution)
print("pcbo solution:", solution)
print("objective function:", H.value(solution), "\n")

print("The solution is", "valid" if H.is_solution_valid(solution) else "invalid")
```

    objective function: 1.0 
    
    quso solution: {0: -1, 1: 1, 2: 1, 3: 1, 4: -1, 5: -1, 6: -1} 
    
    pcbo solution: {'e0': 1, 'e1': 0, 'e2': 0, 'e3': 0, 'z0': 1, 'z1': 1, 'z2': 1}
    objective function: 1 
    
    The solution is valid


Again this matches the result of `H.solve_bruteforce()`.

44acbe7cb1173d0f65a2ab097b06de8b00e61e3e
ipynb
Jupyter Notebook
notebook_examples/FaultTree_example.ipynb
panaali/qubovert
d5ea46349d2a058954fb2cb06f559c0d3fb382c5
Apache-2.0
2020-07-10T20:46:50.000Z
2021-12-29T05:01:55.000Z
notebook_examples/FaultTree_example.ipynb
panaali/qubovert
d5ea46349d2a058954fb2cb06f559c0d3fb382c5
Apache-2.0
2020-06-12T22:37:59.000Z
2022-02-22T20:11:13.000Z
notebook_examples/FaultTree_example.ipynb
panaali/qubovert
d5ea46349d2a058954fb2cb06f559c0d3fb382c5
Apache-2.0
2020-05-13T06:02:38.000Z
2022-03-22T20:45:23.000Z
Qwen/Qwen-72B
1. YES
2. YES