N-Monitoring

Description

The $N$-monitoring problem is described in 1, sections 4.1 and 6.1.

Influence Diagram

The influence diagram of the generalized $N$-monitoring problem where $N\ge 1$ and indices $k=1,...,N$. The nodes are associated with states as follows. Load state $L=\{high,low\}$ denotes the load on a structure, report states $R_k=\{high,low\}$ report the load state to the action states $A_k=\{yes,no\}$ which represent different decisions to fortify the structure. The failure state $F=\{failure,success\}$ represents whether or not the (fortified) structure fails under the load $L$. Finally, the utility at target $T$ depends on the fortification costs and whether F fails.

We begin by choosing $N$ and defining our fortification cost function. We draw the cost of fortification $c_k\sim U(0,1)$ from a uniform distribution, and the magnitude of fortification is directly proportional to the cost. Fortification is defined as

$$\begin{array}{r}f(A_k=yes)=c_k\\ f(A_k=no)=0\end{array}$$

import DecisionProgramming as dp
import numpy as np

dp.activate()

N = 4
np.random.seed(13)
c_k = np.random.random(N)
b = 0.03

def fortification(k, a):
    if not a:
        return c_k[k]
    else:
        return 0

Initialising the influence diagram

We initialise the influence diagram before adding nodes to it.

diagram = dp.InfluenceDiagram()

Adding nodes

Add node $L$ which represents the load on the structure. This node is the root node and thus, has an empty information set. Its states describe the state of the load, they are $high$ and $low$.

L = dp.ChanceNode("L", [], ["high", "low"])
diagram.add_node(L)

The report nodes $R_k$ and action nodes $A_k$ are easily added with a for-loop. The report nodes have node $L$ in their information sets and their states are $high$ and $low$. The actions are made based on these reports, which is represented by the action nodes $A_k$ having the report nodes $R_k$ in their information sets. The action nodes have states $yes$ and $no$, which represents decisions whether to fortify the structure or not.

for i in range(N):
    R = dp.ChanceNode(f"R{i}", ["L"], ["high", "low"])
    diagram.add_node(R)
    A = dp.DecisionNode(f"A{i}", [f"R{i}"], ["yes", "no"])
    diagram.add_node(A)

The failure node $F$ has the load node $L$ and all of the action nodes $A_k$ in its information set. The failure node has states $failure$ and $success$.

F = dp.ChanceNode(
    "F",
    ["L", *[f"A{i}" for i in range(N)]],
    ["failure", "success"]
)
diagram.add_node(F)

The value node $T$ is added as follows.

T = dp.ValueNode("T", ["F", *[f"A{i}" for i in range(N)]])
diagram.add_node(T)

Generating arcs

Now that all of the nodes have been added to the influence diagram we generate the arcs between the nodes. This step automatically orders the nodes, gives them indices and reorganises the information into the appropriate form.

diagram.generate_arcs()

Load State Probabilities

After generating the arcs, the probabilities and utilities can be added. The probability that the load is high, $\mathbb{P}(L=high)$, is drawn from a uniform distribution. For different syntax options for adding probabilities and utilities, see the usage page.

r = np.random.random()
X_L = [r, 1.0-r]
diagram.set_probabilities("L", X_L)

Reporting Probabilities

The probabilities of the report states correspond to the load state. We draw the values $x\sim U(0,1)$ and $y\sim U(0,1)$ from uniform distributions.

$$\begin{array}{r}\mathbb{P}(R_k=high\mid L=high)=max\{x,1-x\}\\ \mathbb{P}(R_k=low\mid L=low)=max\{y,1-y\}\end{array}$$

The probability of a correct report is thus in the range [0.5,1]. (This reflects the fact that a probability under 50% would not even make sense, since we would notice that if the test suggests a high load, the load is more likely to be low, resulting in that a low report “turns into” a high report and vice versa.)

In Decision Programming we add these probabilities by declaring probability matrices for nodes $R_k$. The probability matrix of a report node $R_k$ has dimensions (2,2), where the rows correspond to the states $high$ and $low$ of its predecessor node $L$ and the columns to its own states $high$ and $low$.

for i in range(N):
    x, y = np.random.random(2)
    x = np.max([x, 1-x])
    y = np.max([y, 1-y])
    X_R = diagram.construct_probability_matrix(f"R{i}")
    X_R["high", "high"] = x
    X_R["high", "low"] = 1 - x
    X_R["low", "low"] = y
    X_R["low", "high"] = 1 - y
    diagram.set_probabilities(f"R{i}", X_R)

Probability of Failure

The probability of failure is decreased by fortification actions. We draw the values $x\sim U(0,1)$ and $y\sim U(0,1)$ from uniform distribution.

$$\begin{array}{r}\mathbb{P}(F=failure\mid A_N,\dots ,A_1,L=high)=\frac{max\{x,1-x\}}{exp(b\sum _{k=1,\dots ,N}f(A_k))}\\ \mathbb{P}(F=failure\mid A_N,\dots ,A_1,L=high)=\frac{max\{y,1-y\}}{exp(b\sum _{k=1,\dots ,N}f(A_k))}\end{array}$$

First we initialise the probability matrix for node $F$.

X_F = diagram.construct_probability_matrix("F")

This matrix has dimensions (2, 2, 2, 2, 2, 2) because node $L$ and nodes $A_k$, which form the information set of $F$, all have 2 states and node $F$ itself also has 2 states. The orange colored dimensions correspond to the states of the action nodes $A_k$.

To set the probabilities we have to iterate over the information states. Here it helps to know that in Decision Programming the states of each node are mapped to numbers in the back-end. For instance, the load states $high$ and $low$ are referred to as 1 and 2. The same applies for the action states $yes$ and $no$, they are states 1 and 2. The dp.Diagram.Paths class allows us to iterate over the subpaths of specific nodes. In these paths, the states are referred to by their indices. Using this information, we can easily iterate over the information states using the dp.Diagram.Paths class and enter the probability values into the probability matrix.

x, y = np.random.random(2)
for path in dp.Diagram.Paths([2]*N):
    forticications = [fortification(k, a) for k, a in enumerate(path)]
    denominator = np.exp(b * np.sum(forticications))
    X_F[(0, *path, 0)] = max(x, 1-x) / denominator
    X_F[(0, *path, 1)] = 1.0 - max(x, 1-x) / denominator
    X_F[(1, *path, 0)] = min(y, 1-y) / denominator
    X_F[(1, *path, 1)] = 1.0 - min(y, 1-y) / denominator

After declaring the probability matrix, we add it to the influence diagram.

diagram.set_probabilities("F", X_F)

Utility

The utility from the different scenarios of the failure state at target $T$ are

$$\begin{array}{r}g(F=failure)=0\\ g(F=success)=100.\end{array}$$

Utilities from the action states $A_k$ at target $T$ are

$$\begin{array}{r}g(A_k=yes)=c_k\\ g(A_k=no)=0.\end{array}$$

The total cost is thus

$$Y(F,A_N,\dots ,A_1)=g(F)+(-f(A_N))+\cdots +(-f(A_1))$$

We first declare the utility matrix for node $T$.

Y_T = diagram.construct_utility_matrix('T')

This matrix has dimensions (2, 2, 2, 2, 2) where the dimensions correspond to the numbers of states the nodes in the information set have. Similarly as before, the first dimension corresponds to the states of node $F$ and the other 4 dimensions (in orange) correspond to the states of the $A_k$ nodes. The utilities are set and added similarly to how the probabilities were added above.

for path in dp.Diagram.Paths([2]*N):
    forticications = [fortification(k, a) for k, a in enumerate(path)]
    cost = -sum(forticications)
    Y_T[(0, *path)] = 0 + cost
    Y_T[(1, *path)] = 100 + cost

diagram.set_utility('T', Y_T)

Generate Influence Diagram

The full influence diagram can now be generated. We use the default path probabilities and utilities, which are the default setting in this function. In the Contingent Portfolio Programming example, we show how to use a user-defined custom path utility function.

In this particular problem, some of the path utilities are negative. In this case, we choose to use the positive path utility transformation, which translates the path utilities to positive values. This allows us to exclude the probability cut in the next section.

diagram.generate(positive_path_utility=True)

Decision Model

We initialise the JuMP model and declare the decision and path compatibility variables. Since we applied an affine transformation to the utility function, the probability cut can be excluded from the model formulation.

model = dp.Model()
z = diagram.decision_variables(model)
x_s = diagram.path_compatibility_variables(
    model, z,
    probability_cut=False
)

The expected utility is used as the objective and the problem is solved using Gurobi.

EV = diagram.expected_value(model, x_s)
model.objective(EV, "Max")

model.setup_Gurobi_optimizer(
   ("IntFeasTol", 1e-9),
)
model.optimize()

Analyzing Results

We obtain the decision strategy, state probabilities and utility distribution from the solution.

Z = z.decision_strategy()
S_probabilities = diagram.state_probabilities(Z)
U_distribution = diagram.utility_distribution(Z)

The decision strategy shows us that the optimal strategy is to make all four fortifications regardless of the reports.

In [1]: S_probabilities.print_decision_strategy()

Out[2]:
┌────────────────┬────────────────┐
│ State(s) of R1 │ Decision in A1 │
├────────────────┼────────────────┤
│ high           │ yes            │
│ low            │ yes            │
└────────────────┴────────────────┘
┌────────────────┬────────────────┐
│ State(s) of R2 │ Decision in A2 │
├────────────────┼────────────────┤
│ high           │ yes            │
│ low            │ yes            │
└────────────────┴────────────────┘
┌────────────────┬────────────────┐
│ State(s) of R3 │ Decision in A3 │
├────────────────┼────────────────┤
│ high           │ yes            │
│ low            │ yes            │
└────────────────┴────────────────┘
┌────────────────┬────────────────┐
│ State(s) of R4 │ Decision in A4 │
├────────────────┼────────────────┤
│ high           │ yes            │
│ low            │ yes            │
└────────────────┴────────────────┘

The state probabilities for strategy $Z$ are also obtained. These tell the probability of each state in each node, given strategy $Z$.

In [2]: S_probabilities.print(["L"])

Out[2]:
┌────────┬──────────┬──────────┬─────────────┐
│   Node │     high │      low │ Fixed state │
│ String │  Float64 │  Float64 │      String │
├────────┼──────────┼──────────┼─────────────┤
│      L │ 0.564449 │ 0.435551 │             │
└────────┴──────────┴──────────┴─────────────┘

In [3]: report_nodes = [f"R{i}" for i in range(N)]
In [4]: S_probabilities.print(report_nodes)

Out[4]:
┌────────┬──────────┬──────────┬─────────────┐
│   Node │     high │      low │ Fixed state │
│ String │  Float64 │  Float64 │      String │
├────────┼──────────┼──────────┼─────────────┤
│     R1 │ 0.515575 │ 0.484425 │             │
│     R2 │ 0.442444 │ 0.557556 │             │
│     R3 │ 0.543724 │ 0.456276 │             │
│     R4 │ 0.552515 │ 0.447485 │             │
└────────┴──────────┴──────────┴─────────────┘

In [5]: reinforcement_nodes = [f"A{i}" for i in range(N-1)]
In [6]: S_probabilities.print(reinforcement_nodes)

Out[6]:
┌────────┬──────────┬──────────┬─────────────┐
│   Node │      yes │       no │ Fixed state │
│ String │  Float64 │  Float64 │      String │
├────────┼──────────┼──────────┼─────────────┤
│     A1 │ 1.000000 │ 0.000000 │             │
│     A2 │ 1.000000 │ 0.000000 │             │
│     A3 │ 1.000000 │ 0.000000 │             │
│     A4 │ 1.000000 │ 0.000000 │             │
└────────┴──────────┴──────────┴─────────────┘

In [7]: S_probabilities.print(["F"])

Out[7]:
┌────────┬──────────┬──────────┬─────────────┐
│   Node │  failure │  success │ Fixed state │
│ String │  Float64 │  Float64 │      String │
├────────┼──────────┼──────────┼─────────────┤
│      F │ 0.633125 │ 0.366875 │             │
└────────┴──────────┴──────────┴─────────────┘

We can also print the utility distribution for the optimal strategy and some basic statistics for the distribution.

In [8]: U_distribution.print_distribution()

Out[8]:
┌───────────┬─────────────┐
│   Utility │ Probability │
│   Float64 │     Float64 │
├───────────┼─────────────┤
│ -2.881344 │    0.633125 │
│ 97.118656 │    0.366875 │
└───────────┴─────────────┘

In [8]: U_distribution.print_statistics()

Out[8]:
┌──────────┬────────────┐
│     Name │ Statistics │
│   String │    Float64 │
├──────────┼────────────┤
│     Mean │  33.806192 │
│      Std │  48.195210 │
│ Skewness │   0.552439 │
│ Kurtosis │  -1.694811 │
└──────────┴────────────┘

References

1

Salo, A., Andelmin, J., & Oliveira, F. (2019). Decision Programming for Multi-Stage Optimization under Uncertainty, 1–35. Retrieved from http://arxiv.org/abs/1910.09196