RepairableRBD¶
RepairableRBD ¶
RepairableRBD(edges, components, k=None, input_node=None, output_node=None, on_infeasible_rbd='raise', downtime_cost_rate=0.0, capacity=None, repair_crews=None, maintenance_groups=None, ccf_groups=None)
Bases: RBD
A reliability block diagram of repairable components.
Each component alternates between working and failed: it fails after a
time drawn from its reliability model, and is then repaired, as good
as new, after a time drawn from its repairability model. Components fail
and are repaired independently of one another and of the system state
(unless repair_crews makes them wait for a crew), and a component
keeps running (and can fail) while the system is down. The system is up
whenever its working components connect the input node to the output
node.
Two kinds of analysis are offered:
- Exact long-run (steady-state) metrics, built from each component's
availability
MTTF / (MTTF + MTTR)and failure frequency1 / (MTTF + MTTR):mean_availability,node_availability,system_failure_frequency,mean_up_time,mean_down_time,mean_time_between_failures,expected_cost_rate,capacity_distribution(with node capacities) and the importance measures. A component under age or block replacement enters through its renewal cycle instead (seenode_availability). With fewerrepair_crewsthan components, the long-run values come from the Markov chain of the components and the repair queue, for exponential components (seemean_availability). - Monte-Carlo simulation of a finite window
[0, t_simulation]that starts with every component working:availability(availability over time, criticality measures and, with node capacities, the capacity over time and the fraction of a demand delivered) andcost(the distribution of the window's cost).
Times are in the time unit of the component models, and costs in the currency they are given in.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
edges
|
Iterable[tuple[Hashable, Hashable]]
|
The directed edges of the diagram, e.g.
|
required |
components
|
dict
|
One entry per node other than the input and output nodes, keyed by node name (the input and output nodes never fail: a component given for one is refused, as its edge to the output node, or from the input node, is most likely missing). Each is one of:
|
required |
k
|
dict[Hashable, int]
|
k-out-of-n nodes, as |
None
|
input_node
|
Hashable
|
The input (source) node. By default None: the only node with no incoming edge. If given, it must be that node, or a ValueError is raised. |
None
|
output_node
|
Hashable
|
The output (sink) node. By default None: the only node with no outgoing edge. If given, it must be that node, or a ValueError is raised. |
None
|
on_infeasible_rbd
|
str
|
|
'raise'
|
downtime_cost_rate
|
float
|
Cost per unit time the whole system is down (e.g. lost production), by default 0.0 (not priced). |
0.0
|
capacity
|
dict[Hashable, float or dict]
|
Each node's capacity, keyed by node name, by default None: the
throughput it passes while it is up, a positive number in any unit
(the same for every node), or a dict |
None
|
repair_crews
|
int
|
How many repair crews work on the components, by default None: as
many as are needed, so no work waits. At most this many jobs
proceed at once: a job is what brings a component back up (a repair
or replacement, preventive maintenance that takes time, or a test
that takes time), and the component is down from when the job falls
due. A job that finds every crew busy waits, the component down,
until a crew is free, which then takes the waiting job of the
highest |
None
|
maintenance_groups
|
dict
|
Options for the maintenance groups the components' |
None
|
ccf_groups
|
list of CCFGroup
|
Common-cause groups (see |
None
|
Attributes:
| Name | Type | Description |
|---|---|---|
components |
dict
|
Node name -> the model the node is simulated with: a
|
repairability |
dict
|
Node name -> its time-to-repair model (an |
costs |
dict
|
Node name -> |
downtime_cost_rate |
float
|
The system downtime cost rate. |
capacity |
dict
|
The capacities given, keyed by node name, as floats. |
repair_crews |
int or None
|
The number of repair crews (None: as many as are needed). |
acquisition_costs |
dict
|
Node name -> the one-off cost of buying the unit, for the nodes that
declare a non-zero |
input_node |
Hashable
|
The input node. |
output_node |
Hashable
|
The output node. |
nodes |
list
|
The names of the nodes other than the input and output nodes. |
structure_check |
dict
|
The results of the structural checks made at construction. |
COST_KEYS |
tuple[str, ...]
|
The optional cost keys of a spec dict: |
PER_FAILURE_COST_KEYS |
tuple[str, ...]
|
The cost keys charged per failure, which may be distributions:
|
COMPONENT_SPEC_KEYS |
tuple[str, ...]
|
Every key a spec dict may carry. |
PREVENTIVE_KEYS |
tuple[str, ...]
|
The keys of a |
GROUP_KEYS |
tuple[str, ...]
|
The keys of a maintenance group's options: |
REPAIR_MODELS |
tuple[str, ...]
|
The imperfect repair models a |
INSPECTION_KEYS |
tuple[str, ...]
|
The keys of an |
STANDBY_KEYS |
tuple[str, ...]
|
The keys of a |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a spec dict has an unknown key, or a |
TypeError
|
If a component is not a spec dict, a |
KeyError
|
If a spec dict has no |
Examples:
Two pumps in parallel, each failing on average every 10 hours and taking 1 hour on average to repair:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> pump = {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
>>> pumps = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {"a": pump, "b": pump},
... )
>>> round(pumps.mean_availability(), 4) # 1 - (1 / 11) ** 2
0.9917
>>> result = pumps.availability(t_simulation=50, mc_samples=200, seed=0)
>>> window = result.n_simulations * result.time_simulated_to
>>> round(float(result.system_uptime) / window, 4) # simulated
0.9914
The pair as one node of a larger system, in series with a valve that is replaced instantly at a cost of 250 per failure, while lost production costs 1000 per hour of system downtime:
>>> plant = RepairableRBD(
... [("s", "pumps"), ("pumps", "valve"), ("valve", "t")],
... {
... "pumps": pumps,
... "valve": {
... "reliability": surv.Weibull.from_params([500.0, 2.0]),
... "repairability": "instant",
... "replace_cost": 250.0,
... },
... },
... downtime_cost_rate=1000.0,
... )
>>> round(plant.mean_availability(), 4) # the valve is never down
0.9917
>>> round(plant.expected_cost_rate(), 2) # cost per hour, long run
8.83
A NonRepairable can stand in for a spec dict, and one object can
serve several nodes, since each node gets its own copy:
>>> from repyability import NonRepairable
>>> unit = NonRepairable(
... surv.Exponential.from_params([0.1]),
... surv.Exponential.from_params([1.0]),
... )
>>> same = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {"a": unit, "b": unit},
... )
>>> round(same.mean_availability(), 4)
0.9917
>>> same.components["a"] is same.components["b"]
False
has_costs
property
¶
has_costs
Whether any cost has been declared (per-component or system-wide).
True if some component declares a non-zero "repair_cost",
"replace_cost", "downtime_cost", or preventive-maintenance
or inspection "cost" (a cost distribution always counts), is
on a "calendar" policy (charged its calendar's prices), a
maintenance group has a set-up cost, or downtime_cost_rate is
non-zero: whether running the system costs anything. Costs of 0
price nothing, costs declared inside a nested RepairableRBD do
not count (a nested diagram that carries any is warned of when
this one is built), and neither does an
"acquisition_cost" (a one-off cost, see total_cost). When
nothing is priced, not even an acquisition cost, there is no cost
model to evaluate: the cost analyses (expected_cost_rate,
total_cost, expected_cost and cost) refuse, as
routes.Refused, and the cost of availability's result
is None. With an acquisition cost alone, expected_cost_rate
is 0.0, and cost and that result's cost have no running
cost, exactly, with the acquisition cost beside it.
Returns:
| Type | Description |
|---|---|
bool
|
True if anything is priced. |
acquisition_cost
property
¶
acquisition_cost
The one-off cost of buying the components: the sum of their
"acquisition_cost" (0.0 if none is given). Only this RBD's own
components count, not those inside a nested RepairableRBD.
component_status
property
¶
component_status
Whether each component works, in the simulation in progress.
last_change_planned
property
¶
last_change_planned
Whether the last change of the system's state, in the
simulation in progress, was planned (see next_event).
find_irrelevant_components ¶
find_irrelevant_components()
Return the nodes that cannot affect whether the system works.
A node is irrelevant when it is in no minimal path set, e.g. a node
in parallel with a direct edge, which is a connection that never
fails, or with a path through junctions alone (a
PerfectReliability node, such as a vote point, never fails
either). Whether such a node works never changes whether the
system works, so its Birnbaum and structural importance are zero.
They are found by the reduction to modules, with the junctions
folded in as always working, without listing the path sets. An
irrelevant node is still a node of the RBD (it is in
node_names); a junction is no
component, so it is never listed. The same set is found on
construction and stored in structure_check["irrelevant_nodes"],
with structure_check["has_irrelevant_nodes"].
Returns:
| Type | Description |
|---|---|
set
|
The irrelevant node names; empty if every node is relevant. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If no set of working nodes can reach the output node. |
Examples:
Node "b" is bypassed by the direct edge from "a" to "t":
>>> from repyability import RBD
>>> rbd = RBD([("s", "a"), ("a", "t"), ("a", "b"), ("b", "t")])
>>> rbd.find_irrelevant_components()
{'b'}
get_all_path_sets ¶
get_all_path_sets()
Iterate over every path from the input node to the output node.
A thin wrapper of networkx.all_simple_paths: each path is a list
of node names in order, from the input node to the output node.
These are paths of the graph, not reliability path sets: they ignore
k-out-of-n values (a path through a node with k > 1 does not by
itself make the system work) and need not be minimal. For those, use
get_min_path_sets. The number
of paths can grow exponentially with the size of the diagram, so
avoid exhausting the iterator on very large RBDs.
Returns:
| Type | Description |
|---|---|
Iterator[list[Hashable]]
|
A generator of paths, each a list of node names. |
Examples:
>>> from repyability import RBD
>>> rbd = RBD([("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")])
>>> sorted(rbd.get_all_path_sets())
[['s', 'a', 't'], ['s', 'b', 't']]
get_min_path_sets ¶
get_min_path_sets(include_in_out_nodes=True)
Return the minimal path sets of the RBD.
A path set is a set of nodes whose working is enough for the system
to work; it is minimal when no node can be left out. k-out-of-n
values are accounted for: a minimal path set through a node with
k-out-of-n value k combines path sets reaching k of its
predecessors. The sets are expanded from the diagram's modules (see
system_probability) on first
use and cached; each call returns a new set. Their number can grow
exponentially with redundancy (n stages of duplicated units in
series have 2 ** n), but nothing else needs them: the system
probability, the importance measures and the cut sets are found
without listing them. A junction (a node given
PerfectReliability, such as a k-out-of-n vote point) always
works, so no path set needs it, and none lists it.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
include_in_out_nodes
|
bool
|
Whether each path set includes the input and output nodes, by
default True. Pass False for the components only. (The default
differs from
|
True
|
Returns:
| Type | Description |
|---|---|
set[frozenset[Hashable]]
|
The minimal path sets, each a frozenset of node names. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If there is no path set, i.e. no set of working nodes can reach
the output node (e.g. a |
Examples:
A 2-out-of-3 vote "v" needs two of "a", "b" and "c":
>>> from repyability import RBD
>>> rbd = RBD(
... [
... ("s", "a"), ("s", "b"), ("s", "c"),
... ("a", "v"), ("b", "v"), ("c", "v"), ("v", "t"),
... ],
... k={"v": 2},
... )
>>> paths = rbd.get_min_path_sets(include_in_out_nodes=False)
>>> sorted(sorted(p) for p in paths)
[['a', 'b', 'v'], ['a', 'c', 'v'], ['b', 'c', 'v']]
is_system_working ¶
is_system_working(component_status, method)
Return whether the system works, given which components work.
This is the structure function. It is evaluated through the
diagram's modules (see
system_probability): a
series module works when all of its members do, a parallel one when
any does, a k-out-of-n one when at least k do, and what is left
(e.g. a bridge) through its minimal path sets (method="p": some
path set has every member working) or cut sets (method="c":
every cut set has a member working). Both give the same answer;
"p" is typically faster as it does not need the cut sets. The
input and output nodes need no entry. The layout is worked out once,
so repeated calls (as in the simulations) are cheap.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
component_status
|
dict[Any, bool]
|
Whether each component is working (truthy) or failed (falsy), keyed by node name. Every node in a minimal path set needs an entry; other keys are ignored. |
required |
method
|
str
|
|
required |
Returns:
| Type | Description |
|---|---|
bool
|
True if the system is working, otherwise False. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
KeyError
|
If a node needed for the evaluation has no entry in
|
Examples:
Two parallel nodes "a" and "b" in series with "c":
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> rbd.is_system_working({"a": False, "b": True, "c": True}, "p")
True
>>> rbd.is_system_working({"a": True, "b": True, "c": False}, "c")
False
system_timeline ¶
system_timeline(timelines)
The system's up/down history from its components': the structure function followed over time.
The components' timelines are merged up the diagram's modules (see
system_probability): a
series module is up while all its members are, a parallel one while
any is, a k-out-of-n one while at least k are, and what is left
(e.g. a bridge) while one of its minimal path sets is all up. Each
change of the system's timeline keeps its cause, the component
whose change made it, and whether that change was planned; changes
at the same time are taken in the order of the diagram's
components, each one's in its own order (see
repyability.timelines).
It takes any timelines: an outage log, a what-if edit of one, or
simulated histories (see
RepairableRBD.simulate_timelines).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
timelines
|
Mapping
|
|
required |
Returns:
| Type | Description |
|---|---|
Timeline or Timelines
|
The system's: a Timelines if any component's is one, else a Timeline. Its causes are node names. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a component the system depends on has no timeline, a key is not a component (the input or output node, a repeated node's copy, or an unknown name), or the timelines do not share a window or a number of histories. |
TypeError
|
If |
Examples:
Two pumps in parallel feeding a valve, from their outage logs:
>>> from repyability import RBD, Timeline
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "v"), ("b", "v"), ("v", "t")]
... )
>>> logs = {
... "a": Timeline.from_outages([(10, 30)], end=100),
... "b": Timeline.from_outages([(20, 25), (60, 70)], end=100),
... "v": Timeline.from_outages([(80, 81)], end=100),
... }
>>> plant = rbd.system_timeline(logs)
>>> plant.down_intervals.tolist()
[[20.0, 25.0], [80.0, 81.0]]
>>> plant.failures_by_cause()
{'a': 0, 'b': 1, 'v': 1}
minimal_path_sets ¶
minimal_path_sets()
The minimal path sets of the components (no input or output
node), as a list, smallest first (then by name): the same as
get_min_path_sets
(include_in_out_nodes=False), in the form a
FaultTree's minimal_path_sets
gives, so the same code runs on either.
Returns:
| Type | Description |
|---|---|
list of frozenset
|
The minimal path sets. |
Examples:
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> [sorted(p) for p in rbd.minimal_path_sets()]
[['a', 'c'], ['b', 'c']]
minimal_cut_sets ¶
minimal_cut_sets()
The minimal cut sets of the components, as a list, smallest first
(then by name): the same as
get_min_cut_sets(), in
the form a FaultTree's
minimal_cut_sets gives, so the same code runs on either.
Returns:
| Type | Description |
|---|---|
list of frozenset
|
The minimal cut sets. |
Examples:
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> [sorted(c) for c in rbd.minimal_cut_sets()]
[['c'], ['a', 'b']]
get_min_cut_sets ¶
get_min_cut_sets(include_in_out_nodes=False)
Return the minimal cut sets of the RBD.
A cut set is a set of nodes whose failure is enough for the system
to fail; it is minimal when no node can be left out. The minimal cut
sets are the minimal transversals (hitting sets) of the minimal path
sets, and k-out-of-n values are accounted for. They are built from
the diagram's modules without listing the path sets: a series
module's cut sets are its members', a parallel module's combine one
of each member's, and a k-out-of-n module's combine one of each of
n - k + 1 members'; the part that is not series-parallel (e.g. a
bridge) has its own read off the exact engine's Shannon
decomposition (see minimal_cut_sets_from_path_sets in this
module). They are worked out on first use and cached; each call
returns a new set. A junction (a node given PerfectReliability,
such as a k-out-of-n vote point) never fails, so a cut set with one
never happens, and none is listed.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
include_in_out_nodes
|
bool
|
Whether to derive the cut sets from the path sets including the
input and output nodes, by default False (components only). If
True, the input and output nodes each appear as an extra
single-node cut set. (The default differs from
|
False
|
Returns:
| Type | Description |
|---|---|
set[frozenset[Hashable]]
|
The minimal cut sets, each a frozenset of node names. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If no set of working nodes can reach the output node. |
Examples:
Two parallel nodes "a" and "b" in series with "c":
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> sorted(sorted(c) for c in rbd.get_min_cut_sets())
[['a', 'b'], ['c']]
>>> cuts = rbd.get_min_cut_sets(include_in_out_nodes=True)
>>> sorted(sorted(c) for c in cuts)
[['a', 'b'], ['c'], ['s'], ['t']]
path_set_probabilities ¶
path_set_probabilities(node_probabilities)
Return the probability that each minimal path set fully works.
For each minimal path set (components only) this is the product of
its nodes' probabilities, i.e. the probability that all of them
work, assuming independent nodes. The path sets share nodes, so these
values do not add up to the system probability: use
system_probability for that.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
node_probabilities
|
dict
|
The probability that each node works (not a time: a reliability at some time, or an availability), keyed by node name, as floats or numpy arrays of equal length (not lists). Every node in a minimal path set needs an entry; other keys are ignored. |
required |
Returns:
| Type | Description |
|---|---|
ndarray
|
One value per minimal path set, in the order
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If a node in a minimal path set has no entry. |
Examples:
Two parallel nodes "a" and "b" in series with "c" have
the minimal path sets {a, c} and {b, c}:
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> probs = rbd.path_set_probabilities({"a": 0.9, "b": 0.8, "c": 0.95})
>>> sorted(round(float(p), 4) for p in probs)
[0.76, 0.855]
system_probability ¶
system_probability(node_probabilities, method='p')
Return the exact system probability from each node's probability.
Given the probability that each node works -- e.g. its reliability at some time or its availability -- this returns the probability that the system works, assuming the nodes are independent. It is the engine behind the subclasses' reliability, availability and importance calculations.
The result is exact. On construction the diagram is reduced to
modules, each with a closed form: a series chain works with
probability p1 * p2 * ..., a parallel group with
1 - (1 - p1) * (1 - p2) * ..., and a k-out-of-n group by summing
over how many of its members work. Whatever is not series-parallel
(e.g. a bridge) is left as a core over modules and single nodes, and
is expanded by a Shannon (pivotal) decomposition over its minimal
path sets, with repeated sub-problems solved once. So a
series-parallel diagram is evaluated without listing its path sets,
however many there are (n stages of duplicated units in series
have 2 ** n), and only the core pays the combinatorial price.
Both depend only on the structure, so they are worked out once and
reused.
With method="p" the probability that the system works is
computed, and with method="c" the probability that it fails,
whose complement is returned. Both give the same result (up to
rounding); "p" is the default. Every step is a sum of products
of node probabilities and their complements, so the probability
computed keeps its full relative precision however close to 0 it
is.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
node_probabilities
|
Dict
|
The probability that each node works, keyed by node name: floats, or 1-d arrays (e.g. one value per time) that all have the same length. Every intermediate node needs an entry; entries for the input and output nodes, and any other keys, are not used. The dict is not modified. |
required |
method
|
str
|
|
'p'
|
Returns:
| Type | Description |
|---|---|
ndarray
|
The system probability, a 1-d array with one value per element of the node arrays. It is always an array: float inputs give a one-element array. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
KeyError
|
If an intermediate node has no entry in |
Examples:
Two parallel nodes "a" and "b" in series with "c" work
with probability (1 - 0.1 * 0.2) * 0.95:
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> p = rbd.system_probability({"a": 0.9, "b": 0.8, "c": 0.95})
>>> round(float(p[0]), 4)
0.931
Arrays give one system probability per element:
>>> import numpy as np
>>> p = rbd.system_probability(
... {
... "a": np.array([0.9, 0.5]),
... "b": np.array([0.8, 0.5]),
... "c": np.array([0.95, 0.5]),
... }
... )
>>> [round(float(v), 4) for v in p]
[0.931, 0.375]
system_capacity ¶
system_capacity(node_probabilities)
Return the exact distribution of the system's capacity, from each node's probability of working.
Each node carries its capacity while it works and nothing once
it has failed (a node given none limits nothing), and the system's
capacity is the most that can flow from the input to the output:
the maximum flow of the diagram, which is the least total capacity
of a cut. A series chain carries the least of its nodes' capacities
and a parallel group their sum, and a k-out-of-n node passes flow
only while at least k of its inputs are reached, so the
capacity is positive exactly when the system works. With three
pumps of half the demand each, one failure costs nothing and two
cost half the output.
The distribution is exact, over every combination of the nodes'
states, and assumes they are independent. It is worked out as the
system probability is: each series, parallel or k-out-of-n module's
distribution in closed form from its members' (the least of their
capacities, their sum, their sum while at least k work), and
what is left (e.g. a bridge) by conditioning on its parts'
capacities one at a time. So it costs about what the system
probability does. The subclasses give it the node probabilities:
NonRepairableRBD.capacity_distribution
at a time, and
RepairableRBD.capacity_distribution
in the long run.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
node_probabilities
|
Dict
|
The probability that each node works, keyed by node name: floats, or 1-d arrays (e.g. one value per time) that all have the same length. Every intermediate node needs an entry; other keys are not used. |
required |
Returns:
| Type | Description |
|---|---|
CapacityDistribution
|
The capacities the system can have and their probabilities: one per level for float probabilities, and one row per level and one column per element for arrays. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If no node has a capacity, the diagram is not a valid RBD, the
probability arrays are not all the same length, or a node takes
its capacity from its model (a |
KeyError
|
If an intermediate node has no entry in |
Examples:
Two pumps of 60 and 40 in parallel, feeding a pipe that carries 80:
>>> from repyability import RBD
>>> plant = RBD(
... [("s", "p1"), ("s", "p2"), ("p1", "pipe"), ("p2", "pipe"),
... ("pipe", "t")],
... capacity={"p1": 60, "p2": 40, "pipe": 80},
... )
>>> up = {"p1": 0.9, "p2": 0.8, "pipe": 1.0}
>>> capacity = plant.system_capacity(up)
>>> capacity.levels.tolist()
[0.0, 40.0, 60.0, 80.0]
>>> capacity.probabilities.round(4).tolist()
[0.02, 0.08, 0.18, 0.72]
The pipe, not the pumps, limits the plant when both run.
improvement_allocation ¶
improvement_allocation(target, node_probabilities, fixed=None, weights=None)
Improve the node probabilities just enough to meet a target.
Reliability allocation by a common improvement: every node that is
not fixed has its unreliability 1 - p multiplied by
exp(-x * w), where w is the node's weight and x is one
scale factor shared by all nodes,
p_new = 1 - (1 - p) * exp(-x * w),
and x is solved for (by a bracketed root search) so that
system_probability of the
new probabilities equals target. With the default equal weights
every free node's unreliability shrinks by the same factor, keeping
their ratios. A larger weight changes a node more; a weight of 0
leaves it unchanged, as does a probability of exactly 1.
A target below the current system probability gives x < 0,
which increases the unreliabilities, each capped at 1 (a node
probability of 0): the result is then the lowest node probabilities
that still meet the target. A target the scaling cannot reach
(e.g. because fixed nodes limit the system probability) raises a
ValueError giving the reachable range. A target equal to the best
reachable probability gives the free nodes probability 1. The
solver's result is stored on the RBD as res, a scipy
OptimizeResult whose x holds x (inf when the free
nodes are made perfect), replacing any earlier one.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system probability to reach, in [0, 1]. |
required |
node_probabilities
|
Dict
|
The current probability that each node works, as a single float in [0, 1] per node, keyed by node name. An intermediate node with no entry starts at 0.5. Any other keys (e.g. the input and output nodes) are scaled like the rest and returned. |
required |
fixed
|
list
|
Nodes whose probability is not changed, by default None (every node may change). |
None
|
weights
|
dict
|
A weight per node, by default None (1.0 for every node). If
given, it needs an entry for every node that is not |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
The allocated probability (a float) of every key of
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
KeyError
|
If |
Examples:
Raising two parallel nodes "a" and "b" in series with
"c" from a system probability of 0.931 to 0.99:
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> current = {"a": 0.9, "b": 0.8, "c": 0.95}
>>> new = rbd.improvement_allocation(0.99, current)
>>> {k: round(v, 4) for k, v in sorted(new.items())}
{'a': 0.9814, 'b': 0.9627, 'c': 0.9907}
>>> round(float(rbd.system_probability(new)[0]), 4)
0.99
Every unreliability was cut by the same factor:
>>> sorted({round((1 - new[k]) / (1 - current[k]), 4) for k in new})
[0.1863]
equal_allocation ¶
equal_allocation(target)
Give every node the same probability, chosen to meet a target.
Finds the single probability p that, given to every intermediate
node, makes system_probability
equal target: e.g. target ** (1 / n) for n nodes in
series and 1 - (1 - target) ** (1 / n) for n in parallel.
Any node models are ignored. It runs
improvement_allocation
from 0.5 for every node, with equal weights and nothing fixed, so the
solver's result is stored on the RBD as res. A target of 1 gives
every node probability 1.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system probability to reach, in [0, 1]. |
required |
Returns:
| Type | Description |
|---|---|
dict
|
The allocated probability (a float, the same for every node), keyed by intermediate node name. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
Examples:
Three nodes in series each need 0.9 ** (1 / 3):
>>> from repyability import RBD
>>> rbd = RBD([("s", 1), (1, 2), (2, 3), (3, "t")])
>>> new = rbd.equal_allocation(0.9)
>>> {k: round(v, 4) for k, v in sorted(new.items())}
{1: 0.9655, 2: 0.9655, 3: 0.9655}
simple_allocation ¶
simple_allocation(target, weights=None)
Meet a target with the smallest change in the node log-odds.
Starting from 0.5 for every node, finds the node probabilities that
meet the target with the least weighted change on the log-odds
scale. With s_i = log(p_i / (1 - p_i)) a node's log-odds (0 at
0.5) and w_i its weight, it minimises sum_i s_i ** 2 / w_i
subject to system_probability
equalling target. Any node models are ignored.
At the solution each node's change is proportional to its weight
times the sensitivity of the system's log-odds to it,
p_i (1 - p_i) I_B(i) / (R (1 - R)) with I_B(i) its Birnbaum
importance: nodes that matter more to the system (e.g. a node in
series with a redundant pair), and nodes with larger weights, move
further from 0.5 (at equal sensitivity, twice the weight moves a node
twice as far), and a node with weight 0 stays at 0.5. With equal
weights, nodes placed symmetrically (e.g. all in series, or all in
parallel) get the same value, as
equal_allocation gives.
It is a structural allocation: it uses no component data, only the
diagram, the target and the weights. With every node at 0.5 the
Birnbaum importance equals the
structural_importance, so
a small change moves each node in proportion to its weight times its
structural importance; for larger changes the importances are
re-evaluated at the new probabilities. It is not one of the classic
named methods; see
cost_based_allocation and
minimum_effort_allocation
for those.
It is solved with scipy.optimize.minimize (trust-constr),
from the point that moves every node in proportion to its weight
just far enough to meet the target, with exact gradients; the
system probability is handled on the log-odds scale from both ends
of the exact engine, so it stays exact at any size. A target of 0
or 1 is approached to within 1e-12. The solver's result is stored
on the RBD as res, replacing any earlier one; its x holds
the log-odds of the nodes with a positive weight, in node order.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system probability to reach, in [0, 1]. |
required |
weights
|
dict
|
A non-negative weight per intermediate node, by default None (1.0 for every node). If given, it needs an entry for every intermediate node. |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
The allocated probability (a float) of every intermediate node, keyed by node name. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
KeyError
|
If |
Examples:
With two parallel nodes "a" and "b" in series with "c",
the series node "c" moves furthest:
>>> from repyability import RBD
>>> rbd = RBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")]
... )
>>> new = rbd.simple_allocation(0.99)
>>> {k: round(v, 4) for k, v in sorted(new.items())}
{'a': 0.9395, 'b': 0.9395, 'c': 0.9936}
>>> round(float(rbd.system_probability(new)[0]), 4)
0.99
minimum_effort_allocation ¶
minimum_effort_allocation(target, node_probabilities)
Albert's minimum-effort allocation for a series system.
The minimization-of-effort algorithm (Albert, 1958; MIL-HDBK-338B)
raises a series system to target with the least total effort:
with the current probabilities in ascending order R_1 <= ... <=
R_n, it raises the k least reliable nodes to one common level
R_0 = (target / (R_{k+1} * ... * R_n)) ** (1 / k),
where k is the largest j with R_j below
(target / (R_{j+1} * ... * R_n)) ** (1 / j), and leaves the other
nodes unchanged. The result holds for any effort function the nodes
share that meets Albert's conditions (the effort of raising a
reliability from x to y is non-negative, grows with y
and adds up over successive steps, among others), such as
y - x or log((1 - x) / (1 - y)): the answer does not depend
on which. It needs no solver, so res is not changed.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system probability to reach, in [0, 1]. |
required |
node_probabilities
|
Dict
|
The current probability that each intermediate node works (a single value in [0, 1] per node), keyed by node name. Every intermediate node needs an entry; other keys are not used. |
required |
Returns:
| Type | Description |
|---|---|
dict
|
The allocated probability (a float) of every intermediate node, keyed by node name. A target the system already meets returns the current probabilities. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
References
A. Albert, "A measure of the effort required to increase reliability", Technical Report No. 43, Applied Mathematics and Statistics Laboratory, Stanford University, 1958.
MIL-HDBK-338B, Electronic Reliability Design Handbook, 1998: "Minimization of effort algorithm".
Examples:
The two weakest of four nodes in series are raised together; the others are already good enough:
>>> from repyability import RBD
>>> rbd = RBD([("s", "a"), ("a", "b"), ("b", "c"), ("c", "d"),
... ("d", "t")])
>>> current = {"a": 0.7, "b": 0.8, "c": 0.9, "d": 0.95}
>>> new = rbd.minimum_effort_allocation(0.6, current)
>>> {k: round(v, 4) for k, v in new.items()}
{'a': 0.8377, 'b': 0.8377, 'c': 0.9, 'd': 0.95}
>>> round(float(rbd.system_probability(new)[0]), 4)
0.6
cost_based_allocation ¶
cost_based_allocation(target, node_probabilities, max_probabilities=None, feasibility=None)
Mettas's cost-based allocation: the cheapest way to a target.
Reliability allocation as an optimisation (Mettas, 2000): find the
node probabilities R_i that meet the system target at the least
total cost sum_i c_i(R_i), where each node's cost of improvement
is
c_i(R_i) = exp((1 - f_i) * (R_i - R_min_i) / (R_max_i - R_i)),
R_min_i is its current probability, R_max_i the most it can
reach and f_i in [0, 1) its feasibility: how easily it can be
improved relative to the others (the cost rises faster for a lower
f_i). The cost is 1 at the current probability and grows without
bound towards the maximum, so every node stays within its bounds and
none reaches its maximum. It works on any structure: the system
probability is exact, and the improvement goes where it is cheapest
per unit of system probability, which depends on each node's
importance, current value, maximum and feasibility.
It is solved with scipy.optimize.minimize (SLSQP), from the
point where every node has closed the same fraction of its gap to
its maximum, using exact gradients. The system probability is
handled on the log-odds scale, computed from both ends of the exact
engine, and the cost as the log of its sum, so neither loses
precision nor overflows, however near 0 or 1 the probabilities are.
The allocation returned meets the target. The solver's result is
stored on the RBD as res, replacing any earlier one; its fun
is the log of the total cost of the nodes allowed to change.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system probability to reach, in [0, 1]. |
required |
node_probabilities
|
Dict
|
The current probability that each intermediate node works (a
single value in [0, 1] per node), keyed by node name: the lower
bound |
required |
max_probabilities
|
dict
|
The most each node's probability can be raised to, |
None
|
feasibility
|
dict
|
Each node's feasibility |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
The allocated probability (a float) of every intermediate node, keyed by node name. A target the system already meets returns the current probabilities. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
Warns:
| Type | Description |
|---|---|
UserWarning
|
If the solver stops before converging; the allocation returned still meets the target, but may not be the cheapest. |
References
A. Mettas, "Reliability allocation and optimization for complex systems", Proceedings of the Annual Reliability and Maintainability Symposium, 2000, pp. 216-221.
Examples:
Two parallel pumps in series with a valve, all at 0.9, to reach 0.99. The valve is on every path, so it must reach 0.99 itself; the cheapest allocation takes it only a little further and raises the pumps to supply the rest:
>>> from repyability import RBD
>>> rbd = RBD([("s", "p1"), ("s", "p2"), ("p1", "v"), ("p2", "v"),
... ("v", "t")])
>>> current = {"p1": 0.9, "p2": 0.9, "v": 0.9}
>>> new = rbd.cost_based_allocation(0.99, current)
>>> {k: round(v, 4) for k, v in new.items()}
{'p1': 0.9808, 'p2': 0.9808, 'v': 0.9904}
A valve that is harder to improve stops even closer to 0.99, and the pumps do more:
>>> new = rbd.cost_based_allocation(
... 0.99, current, feasibility={"v": 0.1}
... )
>>> {k: round(v, 4) for k, v in new.items()}
{'p1': 0.9896, 'p2': 0.9896, 'v': 0.9901}
node_names ¶
node_names()
Return the names of the intermediate (component) nodes.
Every node of the diagram except the input and output nodes, in the
order the nodes were added to the graph (for a feasible RBD, the
order of first appearance in edges). Irrelevant nodes are
included. In a NonRepairableRBD a repeated node is the
component it repeats, so only the latter is listed.
Returns:
| Type | Description |
|---|---|
list[Hashable]
|
The intermediate node names, as a new list. |
Examples:
>>> from repyability import RBD
>>> RBD([("s", "a"), ("a", "b"), ("b", "t")]).node_names()
['a', 'b']
structural_importance ¶
structural_importance(working_nodes=None, broken_nodes=None)
Return the structural (probability-free) importance of every node.
The fraction of the states of the other nodes in which the node is
pivotal -- the system works when the node works and fails when it
fails, holding the others fixed. A perfect junction node (a
PerfectReliability drawing device, such as a k-out-of-n vote) is
always working, and is left out. It is computed exactly as the
Birnbaum importance with every node probability at 1/2, so it
depends only on the RBD's structure and not on any failure model --
useful at design time, before any life data exists. It is the same
for a NonRepairableRBD and a
RepairableRBD with the same diagram,
and a common-cause group does not change it.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Iterable[Hashable]
|
Nodes to condition on as working (probability 1), by default None. |
None
|
broken_nodes
|
Iterable[Hashable]
|
Nodes to condition on as failed (probability 0), by default None. With either, the fraction is taken over the states of the nodes that are not forced. A forced node's own importance is still computed, given the other forced nodes. |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float]
|
The structural importance of each intermediate node, in
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If a node is in both |
Examples:
In a two-component parallel system each node is pivotal in half of the other node's states, independent of any failure model:
>>> from surpyval import FixedEventProbability
>>> from repyability import NonRepairableRBD
>>> rbd = NonRepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {
... "a": FixedEventProbability.from_params(0.1),
... "b": FixedEventProbability.from_params(0.4),
... },
... )
>>> si = rbd.structural_importance()
>>> {k: round(v, 4) for k, v in sorted(si.items())}
{'a': 0.5, 'b': 0.5}
Given that "a" has failed, "b" is pivotal in every state
("a"'s own value is unchanged, as it is taken over "b"):
>>> si = rbd.structural_importance(broken_nodes=["a"])
>>> {k: round(v, 4) for k, v in sorted(si.items())}
{'a': 0.5, 'b': 1.0}
to_dict ¶
to_dict()
Serialise the RBD to a JSON-friendly dict.
The dict holds the RBD's type and the repyability_version
that wrote it, plus the constructor inputs as given: the edges, the
node models (or components), the k-out-of-n values and capacities,
the input and output nodes, on_infeasible_rbd, and the
common-cause groups or the downtime cost rate.
from_dict rebuilds the RBD by calling
its constructor again, so the round trip is faithful even for
repeated nodes. Node models are serialised structurally: surpyval's
parametric models in its own format, so an offset, p, f0 and
a fit's covariance round-trip; the RePyability node models (standby,
repeated, load-sharing, regression, NonRepairable,
PerfectReliability and PerfectUnreliability) and nested RBDs
recursively. Per-node
values are stored as lists of entries, so integer and string node
names both survive JSON. Only a
NonRepairableRBD or
RepairableRBD can be serialised.
Returns:
| Type | Description |
|---|---|
dict
|
The serialised RBD. |
Raises:
| Type | Description |
|---|---|
NotImplementedError
|
If a node model cannot be serialised: one that is neither a surpyval model nor a RePyability node model. |
AttributeError
|
If called on a bare |
Examples:
>>> import surpyval as surv
>>> from repyability import NonRepairableRBD
>>> rbd = NonRepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": surv.Weibull.from_params([100, 2])},
... )
>>> d = rbd.to_dict()
>>> d["type"], d["edges"]
('NonRepairableRBD', [['s', 'c'], ['c', 't']])
>>> saved = d["reliabilities"][0]["model"] # surpyval's own format
>>> saved["kind"], saved["model"]["distribution"]
('surpyval', 'Weibull')
to_json ¶
to_json(fp=None, **json_kwargs)
Serialise the RBD to JSON: a string, or written to a file.
Equivalent to json.dumps(self.to_dict(), **json_kwargs); see
to_dict for what is stored. String,
integer and tuple node names all survive: JSON turns a tuple into a
list, and loading turns it back. With fp it is written there, as
surpyval's models' to_json(fp) writes them.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
fp
|
(str, PathLike or file)
|
A path, or a file opened for writing, to write the document to; by default None: it is returned. |
None
|
**json_kwargs
|
Passed to |
{}
|
Returns:
| Type | Description |
|---|---|
str or None
|
The JSON document, or None once written to |
Raises:
| Type | Description |
|---|---|
NotImplementedError
|
If a node model cannot be serialised (see
|
AttributeError
|
If called on a bare |
TypeError
|
If a value cannot be encoded as JSON, e.g. a node name that JSON has no type for. |
Examples:
>>> import json
>>> import surpyval as surv
>>> from repyability import NonRepairableRBD
>>> rbd = NonRepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": surv.Weibull.from_params([100, 2])},
... )
>>> json.loads(rbd.to_json(indent=2))["type"]
'NonRepairableRBD'
from_dict
classmethod
¶
from_dict(d)
Reconstruct an RBD from the output of to_dict.
The RBD is rebuilt by calling its constructor with the stored
inputs, so it is validated again and is equivalent to the original.
Called on a specific subclass, the document's type must match;
called on RBD it dispatches to whichever type the document names.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
d
|
dict
|
A dict made by |
required |
Returns:
| Type | Description |
|---|---|
RBD
|
The reconstructed |
Raises:
| Type | Description |
|---|---|
ValueError
|
If called on a subclass and |
KeyError
|
If a required entry is missing, e.g. |
Examples:
The structure round-trips through a plain dict, so reliability is preserved:
>>> import surpyval as surv
>>> from repyability import NonRepairableRBD
>>> rbd = NonRepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": surv.Weibull.from_params([100, 2])},
... )
>>> restored = NonRepairableRBD.from_dict(rbd.to_dict())
>>> round(restored.sf(50), 4) == round(rbd.sf(50), 4)
True
from_json
classmethod
¶
from_json(s)
Reconstruct an RBD from a JSON document made by to_json.
Equivalent to cls.from_dict(json.loads(s)), so the same type
rules apply (see from_dict). s
can be the document's text or, as surpyval's from_json takes, a
path to the file holding it.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
s
|
(str, PathLike or file)
|
A JSON document made by |
required |
Returns:
| Type | Description |
|---|---|
RBD
|
The reconstructed |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
FileNotFoundError
|
If |
KeyError
|
If a required entry is missing (see |
Examples:
RBD.from_json returns whichever RBD type the document names:
>>> import surpyval as surv
>>> from repyability import RBD, NonRepairableRBD
>>> rbd = NonRepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": surv.Weibull.from_params([100, 2])},
... )
>>> restored = RBD.from_json(rbd.to_json())
>>> type(restored).__name__
'NonRepairableRBD'
>>> round(restored.sf(50), 4)
0.7788
expected_cost_rate ¶
expected_cost_rate(working_nodes=None, broken_nodes=None)
Returns the long-run expected cost per unit time.
Exact (no simulation), from the steady-state quantities:
rate = downtime_cost_rate * (1 - A_sys)
+ sum_i omega_i * (repair_cost_i + replace_cost_i)
+ sum_i nu_i * preventive_cost_i
+ sum_i inspection_cost_i / tau_i
+ sum_i (1 - A_i) * downtime_cost_i
where A_sys is the system's long-run availability
(mean_availability), A_i is node i's (see
node_availability), omega_i is node i's long-run failure
frequency, 1 / (MTTF_i + MTTR_i), and nu_i its frequency of
preventive replacements. So repair_cost and replace_cost
are charged per corrective action (failure), the preventive cost
per preventive action, and the downtime rates per unit time down.
A cost given as a distribution enters through its mean. The result
is in cost per unit of the component models' time.
Under age replacement at T a component renews at a failure or
at a preventive replacement, whichever comes first, in cycles of
mean length C = integral_0^T R + F(T) * MTTR + R(T) * MTTP
(MTTP the mean maintenance time), so omega_i = F(T) / C and
nu_i = R(T) / C: the renewal-reward rate. With instant repair
and maintenance the component's own cost rate is
(c_p * R(T) + c_u * F(T)) / integral_0^T R, as
NonRepairable.cost_rate computes.
Under block replacement at T the renewals are the block times at
which the component is up (see node_availability): with L
the mean time between them and N the mean number of failures in
between, omega_i = N / L and nu_i = 1 / L. With instant
repair and replacement that is (c_p + c_u * M(T)) / T, M the
renewal function of the lives.
A component with hidden failures is inspected every tau_i and,
with a constant failure rate lambda, fails (1 - exp(-lambda *
tau_i)) / tau_i times per unit time (at most once per interval).
With inspection cost c_i and downtime cost rate c_d it
costs c_i / tau + c_d * U(tau), U its unavailability (see
node_availability): frequent tests cost more, and rare ones
leave failures hidden for longer. The rate is least near tau =
sqrt(2 * c_i / (lambda * c_d)).
With limited repair_crews, A_sys, A_i and omega_i
come from the crews' Markov chain (see mean_availability): a
component waiting for a crew is down, and pays its downtime cost,
and it fails at its constant rate while it is up.
Every cost is optional and defaults to 0, so any subset can be
priced; with nothing priced, not even an acquisition cost (see
has_costs), there is no cost model to evaluate, and this
refuses (with an acquisition cost alone it is 0.0). Costs
are undiscounted, and only this RBD's own costs count: costs
declared inside a nested RepairableRBD are left out (a nested
diagram that carries any is warned of when this one is built).
The simulated counterpart is cost, whose cost_rate
converges to this value as the window grows.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working: they never fail, so they incur no corrective or downtime cost, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed: they incur no corrective cost (they never change state) but are down for all time, by default None. |
None
|
Returns:
| Type | Description |
|---|---|
float
|
Expected cost per unit time. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If nothing is priced, not even an acquisition cost (see
|
NotImplementedError
|
If a component is under block replacement with models its
exact values do not cover (see |
Examples:
One component with MTTF 10 and MTTR 1 is up 10/11 of the time and
fails 1/11 times per unit time, so at 100 per repair and 50 per unit
time of outage the cost rate is (100 + 50) / 11:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... "repair_cost": 100.0,
... }
... },
... downtime_cost_rate=50.0,
... )
>>> round(rbd.expected_cost_rate(), 4)
13.6364
total_cost ¶
total_cost(horizon, working_nodes=None, broken_nodes=None, *, discount_rate=0.0)
Returns the total cost of owning the system for horizon.
The life-cycle cost: buying the components, then running the system
for horizon at the long-run cost rate,
total = acquisition_cost + expected_cost_rate() * horizon
with acquisition_cost the sum of the components'
"acquisition_cost". With a discount_rate r it is the
present value: the components are bought at the start, and
the running costs, spent at a steady rate, are discounted
continuously, so the horizon counts as (1 - exp(-r * horizon)) /
r. The running cost is the long-run rate,
exact over a horizon long compared with the components' cycles;
expected_cost(horizon).total is the exact expected cost of
owning the system from new, and cost simulates a window from new
(whose CostResult gives the same acquisition_cost
separately).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
horizon
|
float or array - like
|
How long the system is owned, in the time unit of the component
models: non-negative, and finite unless the costs are discounted
( |
required |
working_nodes
|
Collection[Hashable]
|
As for |
None
|
broken_nodes
|
Collection[Hashable]
|
As for |
None
|
discount_rate
|
float
|
The continuous discount rate per unit time of the component
models, by default 0 (undiscounted): for an annual rate of 7%
with models in hours, |
0.0
|
Returns:
| Type | Description |
|---|---|
float or ndarray
|
The total cost over the horizon, or over each. |
Raises:
| Type | Description |
|---|---|
ValueError
|
As for |
NotImplementedError
|
As for |
Examples:
A pump bought for 20,000, failing on average every 1000 hours, repaired in 10 at 500 per repair, with lost production at 100 per hour, owned for ten years:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "p"), ("p", "t")],
... {
... "p": {
... "reliability": surv.Exponential.from_params([1e-3]),
... "repairability": surv.Exponential.from_params([0.1]),
... "repair_cost": 500.0,
... "acquisition_cost": 20000.0,
... }
... },
... downtime_cost_rate=100.0,
... )
>>> round(rbd.expected_cost_rate(), 4) # (500 + 100 * 10) / 1010
1.4851
>>> round(rbd.total_cost(87600.0))
150099
Its present value at 7% a year, the ten years counting as 63,656 hours:
>>> import math
>>> r = math.log(1.07) / 8760 # per hour
>>> round(rbd.total_cost(87600.0, discount_rate=r))
114538
spares_demand ¶
spares_demand(horizon, *, nodes=None, fleet=1, method='exact', mc_samples=None, seed=None, parts=None)
How many spares each component uses over [0, horizon), from
new: the distribution of its replacements, for one system or a
fleet of fleet; and how many each part that several components
use (parts) takes from its one shelf.
A component uses a spare at each failure and at each preventive
replacement; a standby group at each of its units' failures. Its
replacements are a renewal process: each up time (its life, or the
replacement age if that comes first) ends in one, then a repair or
maintenance time, and it starts again as new. So the n-th
replacement comes after n - 1 whole cycles and an up time, and
the probability of n or more by the horizon is that of their
sum falling before it (a replacement at the horizon itself falls
after it, as in the simulation). method="exact" works that out
on a grid, to about 1e-6, for components with corrective repair
alone or under age replacement. Under block replacement, a unit's
next replacement is at its failure or at the next block time,
whichever comes first (a unit down then skips it), counted block
interval by block interval on the grid. With hidden failures, the
replacements fall on the tests that find them, and are counted
there: exactly, tested and repaired in no time by tests that find
every failure, and otherwise from the cycle between those tests,
followed test by test (with tests that can miss a failure, over
where each cycle starts between the full tests). A fleet's count
is the sum of its systems', independent and each from new.
method="simulate" counts them in mc_samples simulations of
the whole system instead, which also covers standby groups, repair
crews and tests that can last as long as their interval. With
constant failure rates and instant repair, the counts are Poisson.
Identical parts often share one shelf: the seals of a station's
three pumps come from one bin. parts={part: [nodes]} pools the
spares of each part's components: their counts are
independent, so the part's is their sum (by simulation, the sum in
each simulation). Its components may differ, one under age
replacement and the others not.
The exact count convolves a component's replacement times on a grid, once for each replacement likely in the horizon, so it grows with the horizon over the component's life (a few seconds for some hundreds of replacements), and under block replacement with the number of block intervals too: count a long horizon by simulation.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
horizon
|
float
|
The time counted over, from new. |
required |
nodes
|
Collection[Hashable]
|
The components to count for, by default every component, or
none when |
None
|
fleet
|
int
|
How many systems, by default 1. |
1
|
method
|
str
|
|
'exact'
|
mc_samples
|
int
|
The number of simulations for |
None
|
seed
|
int
|
Seeds the simulations, by default None. |
None
|
parts
|
dict
|
|
None
|
Returns:
| Type | Description |
|---|---|
dict[Hashable, SparesDemand]
|
Per component, and per part, the distribution of the spares it
uses, with its |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
For a part with two members in one common-cause group; and
with |
Examples:
A pump with a constant failure rate of 0.01 an hour, replaced in no time, for a fleet of 5 over 1,000 hours: Poisson with mean 50.
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "pump"), ("pump", "t")],
... {
... "pump": {
... "reliability": surv.Exponential.from_params([0.01]),
... "repairability": "instant",
... }
... },
... )
>>> demand = rbd.spares_demand(1000.0, fleet=5)["pump"]
>>> round(demand.mean, 4)
50.0
>>> demand.stock(0.95)
62
spares_stock ¶
spares_stock(lead_time, *, fill_rate=None, stockout_probability=None, nodes=None, fleet=1, parts=None)
The fewest spares of each component (or each part several
share, parts) to hold for a fill rate, or to run out no more
than a fraction of the time, when each spare used is reordered at
once and arrives lead_time later (one-for-one, or (S - 1,
S), replenishment), in the long run.
With a stock of S a spare is on the shelf while fewer than
S are on order: those used in the last lead time. The demand is
the component's replacements (see spares_demand), in the long
run: from a random time the first comes at the end of what is left
of an up time, or after what is left of a repair or maintenance and
an up time, and the next after whole cycles. So the stock-out
probability is that of S or more replacements in a lead time
from a random time, and the fill rate that of fewer than S in
the lead time before a replacement. For a fleet, the systems'
demands add up. Worked out on a grid, to about 1e-6, for components
with corrective repair alone or under age replacement; and on its
tests, for a component with hidden failures (from a random time,
the next replacement is j tests on with probability P(C >=
j) / S, C a cycle's length in tests and S its mean):
exactly, tested and repaired in no time by tests that find every
failure, and otherwise from the cycle followed test by test (with
tests that can miss a failure, over where each cycle starts between
the full tests, see _spares).
Under block replacement every T, with repairs and block
replacements in no time, each block interval starts with a new
unit, so the demand repeats every interval, and a lead time's
depends on where in the interval it starts: it is averaged over
that, uniform on [0, T) from a random time, and as the
replacements fall from a replacement (a failure, at the renewal
density, or a block replacement), on a grid of the interval, to
about 1e-6. With repairs or block replacements that take time, a
unit down at a block time is not replaced there, and an interval
need not start new: the demand is counted from a typical
replacement in the long run, a failure at each phase of the
interval or a block replacement, each followed block interval by
block interval, the replacements before one as those after it,
and from a random time by Campbell's formula; on grids of the
interval, extrapolated, to about 1e-6. A unit dead on arrival
fails as it is put in: with repairs and block replacements in no
time, each replacement brings a geometric number more at its
instant, which the lasting units' count is compounded with
(#322); taking time, they are followed on the grid. A fleet's
systems are taken as on block schedules of their own, out of step
with each other.
A part's components draw on one shelf, which needs fewer
spares than a shelf each: its spares on order are the sum of its
components' independent ones, and a demand comes from component
i with the share rate_i / sum(rates) of their long-run
replacement rates, finding i's as its own demands do and the
others' as at a random time (see the spares
guide).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
lead_time
|
float
|
The time a spare takes to arrive once ordered. |
required |
fill_rate
|
float
|
The fraction of demands to meet from the shelf, in |
None
|
stockout_probability
|
float
|
The most fraction of time to be out of stock, in |
None
|
nodes
|
Collection[Hashable]
|
The components, by default every component (none when
|
None
|
fleet
|
int
|
How many systems share the stock, by default 1. |
1
|
parts
|
dict
|
|
None
|
Returns:
| Type | Description |
|---|---|
dict[Hashable, SparesStock]
|
Per component, and per part, the stock and what it achieves, with the distributions behind them. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
For a part with two members in one common-cause group, or two under block replacement, a component under block replacement whose repairs or block replacements take time or whose life may end at 0, with hidden failures whose tests can last as long as their interval, a standby group or a life that may never end, or while a component can wait for a repair crew, or if more than 2,000 replacements are likely in a lead time. |
Examples:
A pump with a constant failure rate of 0.01 an hour, replaced in no time, and 300 hours to get a spare: 3 are on order on average, and 7 on the shelf meet 96.6% of demands.
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "pump"), ("pump", "t")],
... {
... "pump": {
... "reliability": surv.Exponential.from_params([0.01]),
... "repairability": "instant",
... }
... },
... )
>>> stock = rbd.spares_stock(300.0, fill_rate=0.95)["pump"]
>>> stock.stock
7
>>> round(stock.fill_rate, 4)
0.9665
allocate_redundancy ¶
allocate_redundancy(horizon, *, nodes=None, trains=None, min_availability=None, max_units=None, method='exact', discount_rate=0.0)
Choose the redundancy with the lowest total cost of ownership.
How many identical copies of each node to fit in active parallel
(or of each train of nodes to add alongside it, see trains) so
that owning the system for horizon costs least: each copy costs
its "acquisition_cost" to buy and its running costs (repairs,
replacements, preventive maintenance, inspections, its own downtime
cost) to keep, and together the copies save the cost of the system
being down (downtime_cost_rate). That is total_cost for the
system with the copies drawn out:
total = sum_i n_i * (a_i + horizon * r_i)
+ horizon * downtime_cost_rate * (1 - A_sys)
+ the cost of the nodes not considered
with n_i copies of node i, each bought for a_i and
running at r_i per unit time, and A_sys the system's
long-run availability. With a discount_rate r it is the
present value, as total_cost gives it: the copies are bought at
the start and the running costs discounted, horizon counting as
(1 - exp(-r * horizon)) / r. The copies fail and are repaired
independently, so n copies of a node of long-run availability
A are all down a fraction (1 - A) ** n of the time, and each
design is scored exactly, with no simulation. (Copies of a
component with hidden failures are inspected together, and are
scored over the inspection period, as mean_availability does.)
More copies cost more and save less and less downtime, so the total
is not monotone in them. method="exact" finds a proven optimum
from the greedy solution: a design no worse than that spends no more
on copies than its total, which caps every node's copies, and
(without min_availability) no node takes a copy that could not
pay for itself, since the k + 1-th copy of a node of
unavailability U can save at most horizon *
downtime_cost_rate * U ** k * (1 - U). When every node considered
is in series with the rest of the system (it lies on every path)
and has no hidden failures, the system's availability is theirs
times the rest's, and a dynamic program over the nodes finds the
optimum however many there are; otherwise a branch and bound over
the designs does, which suits a handful of nodes.
method="greedy" adds or removes one copy at a time while that
lowers the total: fast, but not guaranteed optimal.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
horizon
|
float
|
How long the system is owned, in the time unit of the component models: finite and non-negative. |
required |
nodes
|
Collection[Hashable]
|
The components that may be given copies, by default every one
with an |
None
|
trains
|
dict
|
Trains of components that may be given copies, |
None
|
min_availability
|
float
|
Only consider designs whose long-run availability is at least this, in (0, 1), by default no limit. |
None
|
max_units
|
int or dict
|
The most copies (at least 1) of every node or train considered (an int) or of particular ones (a dict; those it leaves out are unlimited), by default unlimited. A node or train whose copies cost nothing over the horizon needs one. |
None
|
method
|
str
|
|
'exact'
|
discount_rate
|
float
|
The continuous discount rate per unit time of the component
models, by default 0 (undiscounted), as for |
0.0
|
Returns:
| Type | Description |
|---|---|
TotalCostAllocation
|
The chosen |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
If a component is under block replacement with models its exact
values do not cover (see |
Examples:
A pump that fails on average every 1000 hours and takes 10 to repair, bought for 20,000 and repaired for 500, when an hour without pumping costs 100. Over ten years (87,600 hours) a second pump pays for itself; a third would not:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "pump"), ("pump", "t")],
... {
... "pump": {
... "reliability": surv.Exponential.from_params([1e-3]),
... "repairability": surv.Exponential.from_params([0.1]),
... "repair_cost": 500.0,
... "acquisition_cost": 20000.0,
... }
... },
... downtime_cost_rate=100.0,
... )
>>> round(rbd.total_cost(87600.0)) # one pump
150099
>>> best = rbd.allocate_redundancy(87600.0)
>>> best.units, round(best.total_cost)
({'pump': 2}, 127591)
Over one year the second pump does not pay for itself:
>>> rbd.allocate_redundancy(8760.0).units
{'pump': 1}
A station needing two of three pump trains, each a pump and its motor, when an hour without pumping costs 2,000: is a fourth train worth it? Name one train, and its copies join the vote:
>>> pump = {
... "reliability": surv.Exponential.from_params([1e-3]),
... "repairability": surv.Exponential.from_params([0.1]),
... "repair_cost": 500.0,
... "acquisition_cost": 20000.0,
... }
>>> motor = {
... "reliability": surv.Exponential.from_params([2e-4]),
... "repairability": surv.Exponential.from_params([0.05]),
... "acquisition_cost": 10000.0,
... }
>>> trains = (1, 2, 3)
>>> station = RepairableRBD(
... [("s", f"pump {i}") for i in trains]
... + [(f"pump {i}", f"motor {i}") for i in trains]
... + [(f"motor {i}", "t") for i in trains],
... {
... **{f"pump {i}": pump for i in trains},
... **{f"motor {i}": motor for i in trains},
... },
... k={"t": 2},
... downtime_cost_rate=2000.0,
... )
>>> round(station.total_cost(87600.0)) # three trains
319927
>>> best = station.allocate_redundancy(
... 87600.0, trains={"train 1": ["pump 1", "motor 1"]}
... )
>>> best.units, round(best.total_cost)
({'train 1': 2}, 295306)
availability_allocation ¶
availability_allocation(target, method='cost_based', *, fixed=None, weights=None, max_availabilities=None, feasibility=None)
What availability each component needs for the system to meet a target, and the MTTF or MTTR that gives it.
Availability allocation: one of the reliability allocation methods
(see cost_based_allocation), run on the components' long-run
availabilities (node_availability), with the system scored as
mean_availability scores it. A component's availability is
MTTF / (MTTF + MTTR), so the availability A allocated to it
is met by an MTTF of MTTR * A / (1 - A) at its current MTTR, or
by an MTTR of MTTF * (1 - A) / A at its current MTTF (or by any
pair in the same ratio): both are reported. To choose the cheapest
pair instead, use mttf_mttr_allocation.
Only components with corrective repair alone are allocated an
availability: those that fail and take time to repair, with no
"preventive" or "inspection" schedule. The others keep
theirs: components with a schedule (whose intervals
optimal_replacement_intervals and
optimal_inspection_intervals choose), components repaired
instantly (never down), nested RepairableRBD nodes, components
with units that never fail or repairs that may never end, the
members of common-cause groups (whose MTTF and MTTR are their
group's; the system is scored over the groups' joint states), and
the components in fixed. A held component that is inspected or
block-replaced is up with a probability that varies over its
schedule, together with the others on the same calendar; it enters
with that variation, as it does in mean_availability, so that
the allocation meets the target exactly.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system's long-run availability to reach, in [0, 1]. |
required |
method
|
str
|
How to allocate it, from the current availabilities:
|
'cost_based'
|
fixed
|
Collection[Hashable]
|
Components that keep their availability. |
None
|
weights
|
dict
|
With |
None
|
max_availabilities
|
dict
|
With |
None
|
feasibility
|
dict
|
With |
None
|
Returns:
| Type | Description |
|---|---|
AvailabilityAllocation
|
The |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
KeyError
|
If |
NotImplementedError
|
As for |
Examples:
Two pumps in parallel, each with an MTTF of 10 h and an MTTR of 1 h, in series with a valve with an MTTF of 50 h and an MTTR of 2 h, are up 95.4% of the time. For 98%, the cheapest allocation asks most of the valve, which every path goes through:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def unit(mttf, mttr):
... return {
... "reliability": surv.Exponential.from_params([1 / mttf]),
... "repairability": surv.Exponential.from_params([1 / mttr]),
... }
>>> plant = RepairableRBD(
... [("s", "p1"), ("s", "p2"), ("p1", "v"), ("p2", "v"),
... ("v", "t")],
... {"p1": unit(10, 1), "p2": unit(10, 1), "v": unit(50, 2)},
... )
>>> round(plant.mean_availability(), 4)
0.9536
>>> need = plant.availability_allocation(0.98)
>>> {node: round(a, 4) for node, a in need.availability.items()}
{'p1': 0.9318, 'p2': 0.9318, 'v': 0.9846}
The valve needs an MTTF of 128 h at its 2 h repairs, or repairs of 0.78 h at its 50 h MTTF:
>>> round(need.mttf["v"]), round(need.mttr["v"], 2)
(128, 0.78)
mttf_mttr_allocation ¶
mttf_mttr_allocation(target, *, levers='both', fixed=None, max_mttf=None, min_mttr=None, mttf_feasibility=None, mttr_feasibility=None)
The cheapest MTTFs and MTTRs that meet a system availability target.
A component's availability, MTTF / (MTTF + MTTR), can be raised
with a longer MTTF (reliability) or a shorter MTTR
(maintainability), which cost different amounts. This is Mettas's
cost-based allocation (see cost_based_allocation) with both
levers: it finds each component's MTTF and MTTR that meet the target
at the least total cost, where lowering its failure rate
lambda = 1 / MTTF from lambda_0 towards its least,
lambda_min = 1 / max_mttf, costs
exp((1 - f) * (lambda_0 - lambda) / (lambda - lambda_min)),
and cutting its MTTR from MTTR_0 towards its least,
min_mttr, costs
exp((1 - f) * (MTTR_0 - MTTR) / (MTTR - min_mttr)),
each with its own feasibility f in [0, 1): the lower it is, the
faster the cost rises. Each cost is 1 while its lever is unused and
grows without bound towards its limit. Without limits, raising the
MTTF by a factor k costs exp((1 - f) * (k - 1)), and so
does cutting the MTTR by one: the availability depends on the ratio
of the two alone, so at equal feasibilities both are used alike.
The cheaper lever is used first, and the dearer one only once the
cheaper one costs as much at the margin.
With levers="mttr" the failure behaviour is held and only the
repairs are allocated (maintainability allocation): a shorter
repair is worth most on a component that is often down, so, other
things equal, the components that fail most often get the shortest
repairs. With levers="mttf" only the MTTFs change.
The components allocated, and those that keep their availability,
are as for availability_allocation, and the system is scored as
mean_availability scores it. It is solved with
scipy.optimize.minimize (SLSQP) with exact gradients, from the
point where every lever has closed the same fraction of its gap to
its limit, and the design returned meets the target. The solver's
result is stored on the RBD as res, replacing any earlier one;
its fun is the log of the total cost of the levers that may
change.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
target
|
float
|
The system's long-run availability to reach, in [0, 1]. |
required |
levers
|
str
|
|
'both'
|
fixed
|
Collection[Hashable]
|
Components that keep their MTTF and MTTR. |
None
|
max_mttf
|
dict
|
The most a component's MTTF can reach, at least its current MTTF (which holds it). A component without one has no limit. |
None
|
min_mttr
|
dict
|
The least a component's MTTR can reach, from 0 (the default for a component without one) to its current MTTR (which holds it). |
None
|
mttf_feasibility
|
dict
|
The feasibility of raising a component's MTTF, in [0, 1), 0.5 without one. |
None
|
mttr_feasibility
|
dict
|
The feasibility of cutting a component's MTTR, in [0, 1), 0.5 without one. |
None
|
Returns:
| Type | Description |
|---|---|
AvailabilityAllocation
|
Every component's |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
As for |
Warns:
| Type | Description |
|---|---|
UserWarning
|
If the solver stops before converging; the design returned still meets the target, but may not be the cheapest. |
Examples:
The pumps and valve of availability_allocation, to 98%: at
equal feasibilities each component's MTTF rises by the factor its
MTTR falls by, and the valve's the most:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def unit(mttf, mttr):
... return {
... "reliability": surv.Exponential.from_params([1 / mttf]),
... "repairability": surv.Exponential.from_params([1 / mttr]),
... }
>>> plant = RepairableRBD(
... [("s", "p1"), ("s", "p2"), ("p1", "v"), ("p2", "v"),
... ("v", "t")],
... {"p1": unit(10, 1), "p2": unit(10, 1), "v": unit(50, 2)},
... )
>>> design = plant.mttf_mttr_allocation(0.98)
>>> {node: round(t, 1) for node, t in design.mttf.items()}
{'p1': 10.6, 'p2': 10.6, 'v': 85.5}
>>> {node: round(t, 2) for node, t in design.mttr.items()}
{'p1': 0.94, 'p2': 0.94, 'v': 1.17}
>>> round(design.system_availability, 4)
0.98
Holding the failure behaviour, only the repairs change:
>>> repairs = plant.mttf_mttr_allocation(0.98, levers="mttr")
>>> {node: round(t, 2) for node, t in repairs.mttr.items()}
{'p1': 0.74, 'p2': 0.74, 'v': 0.78}
optimal_replacement_intervals ¶
optimal_replacement_intervals(nodes=None, *, allowed=None, min_availability=None, max_cost_rate=None, assume_unlimited_crews=False)
Choose the age-replacement intervals for the system as a whole.
NonRepairable.find_optimal_replacement chooses one unit's
replacement age on its own. In a system the components should be
chosen together: a unit whose failure stops the system is worth
replacing sooner than one with a standby, and replacing a unit takes
the system down if its replacement takes time and nothing covers
for it. This chooses the age-replacement interval of each component
in nodes for:
- the lowest long-run cost rate (
expected_cost_rate), by default; - the lowest cost rate that keeps the system's long-run availability
(
mean_availability) at leastmin_availability; - or the highest availability within a cost rate of
max_cost_rate.
The long-run values are exact (see expected_cost_rate), so the
choice is too, to the precision of the search: a gradient search
(SLSQP) over the logarithms of the intervals, from several starting
points, keeping the best. Each interval ranges from a thousandth to a
thousand times the component's mean life; one found at the top of
that range is compared with never replacing the component (an
interval of inf), which is taken if no worse. The cost rate is
usually flat near its minimum, so intervals some way from the ones
found cost almost the same.
Replacements are often made on a calendar (quarterly, yearly): the
intervals are then chosen from allowed, as
optimal_inspection_intervals chooses tests', every combination
tried when there are at most 2000 of them, which gives the
optimum, and a local search (one interval changed at a time, from
several starting points) made otherwise.
The plan is the best for the long run. A plant starts with every
unit new, so the units' first replacements fall due together, and
where they take the system down together (redundant units whose
replacements take time) its first years cost more than the plan's
cost rate says: expected_cost gives the cost from new, and from
units of other ages (state) a staggered start.
With limited repair_crews a component can wait for a crew, and
the exact long-run values do not hold: the choice is refused unless
assume_unlimited_crews, which chooses the intervals as if
every repair started at once. Simulate that plan with the crews,
with_intervals(plan).cost(), to see what the waiting costs.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
nodes
|
Collection[Hashable]
|
The components whose intervals to choose, each under age
replacement (a |
None
|
allowed
|
sequence of float or dict
|
The intervals to choose from: one sequence for every node, or a
dict of a sequence per node; |
None
|
min_availability
|
float
|
The least long-run system availability allowed, in (0, 1). |
None
|
max_cost_rate
|
float
|
The highest long-run cost rate allowed: the intervals then give the highest availability within it. |
None
|
assume_unlimited_crews
|
bool
|
With limited |
False
|
Returns:
| Type | Description |
|---|---|
MaintenancePlan
|
The interval of each component in |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a node is not a component under age replacement; if nothing
is priced; if both targets are given, or one is out of range;
if |
NotImplementedError
|
If the long-run values are not known exactly (see
|
Examples:
A pump that wears out, in series with a pair of them in parallel; repairs take about 23 hours and replacements about 7, and the plant loses 500 an hour while it is down:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def pump():
... return {
... "reliability": surv.Weibull.from_params([1000, 2.5]),
... "repairability": surv.LogNormal.from_params([3.0, 0.5]),
... "replace_cost": 5000.0,
... "preventive": {
... "interval": 1000.0,
... "duration": surv.Weibull.from_params([8, 3]),
... "cost": 1000.0,
... },
... }
>>> plant = RepairableRBD(
... [("s", "a"), ("a", "b1"), ("a", "b2"), ("b1", "t"),
... ("b2", "t")],
... {"a": pump(), "b1": pump(), "b2": pump()},
... downtime_cost_rate=500.0,
... )
>>> plan = plant.optimal_replacement_intervals()
>>> {node: round(t) for node, t in plan.intervals.items()}
{'a': 590, 'b1': 497, 'b2': 497}
>>> round(plan.cost_rate, 2), round(plan.availability, 4)
(20.12, 0.9803)
The pump alone in the line is replaced later than those with a standby: its replacements stop the plant too.
optimal_inspection_intervals ¶
optimal_inspection_intervals(nodes=None, *, allowed=None, min_availability=None, max_cost_rate=None, assume_unlimited_crews=False, offset_shares=None)
Choose the proof-test intervals of components with hidden failures, for the system as a whole.
Testing a component with hidden failures more often costs more
tests, and testing it less often leaves its failures hidden for
longer (see node_availability). This chooses the inspection
interval of each component in nodes for:
- the lowest long-run cost rate (
expected_cost_rate: tests, repairs and downtime), by default; - the lowest cost rate that keeps the system's long-run availability
at least
min_availability: for a safety function, whose unavailability is its average probability of failure on demand, a PFDavg of at most1 - min_availability; - or the highest availability within a cost rate of
max_cost_rate.
Tests are usually made on a calendar (monthly, quarterly, yearly),
and components tested at the same times are down together, so the
long-run values depend on how the schedules line up. So the
intervals are chosen from allowed: every combination is tried
when there are at most 2000 of them, which gives the optimum, and a
local search (changing one interval at a time, from several
starting points) is made otherwise. With allowed left out, the
RBD must have one component with hidden failures, whose interval
is then searched continuously, as by
optimal_replacement_intervals.
When redundant components are tested matters as well: two tested
at once are down together for as long as a failure of both stays
hidden, where tests half an interval apart find a common-cause
failure twice as soon. offset_shares chooses the times of the
first tests with the intervals.
A component whose tests can miss a failure (a coverage below
1) keeps its full tests' interval (full_test), a whole number
of its test intervals: its interval is chosen among those that
divide it, from allowed (which must) or, left out, from every
one in range.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
nodes
|
Collection[Hashable]
|
The components whose intervals to choose, each with hidden
failures (an |
None
|
allowed
|
sequence of float or dict
|
The intervals to choose from: one sequence for every node, or a dict of a sequence per node. Left out, the one component with hidden failures has its interval chosen from all. |
None
|
min_availability
|
float
|
The least long-run system availability allowed, in (0, 1). |
None
|
max_cost_rate
|
float
|
The highest long-run cost rate allowed: the intervals then give the highest availability within it. |
None
|
assume_unlimited_crews
|
bool
|
With limited |
False
|
offset_shares
|
sequence of float, dict or str
|
Choose each node's offset, the time of its first test, as well:
as a share of its interval, in [0, 1), from these, one sequence
for every node or a dict of one per node. |
None
|
Returns:
| Type | Description |
|---|---|
MaintenancePlan
|
The interval of each component in |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a node is not a component, or has no hidden failures; if
|
NotImplementedError
|
If a component's hidden failures have no exact long-run values
(see |
Examples:
A shutdown valve whose dangerous failures are hidden, at 2e-6
per hour, each proof test costing 500: the cheapest monthly,
quarterly, half-yearly, yearly or two-yearly tests (in hours) that
keep the PFDavg at most 1e-3, for one valve and for two in
parallel (1oo2):
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def valve():
... return {
... "reliability": surv.Exponential.from_params([2e-6]),
... "repairability": "instant",
... "inspection": {"interval": 8760.0, "cost": 500.0},
... }
>>> calendar = [730.0, 2190.0, 4380.0, 8760.0, 17520.0]
>>> one = RepairableRBD([("s", "v"), ("v", "t")], {"v": valve()})
>>> one.optimal_inspection_intervals(
... allowed=calendar, min_availability=1 - 1e-3
... ).intervals
{'v': 730.0}
>>> pair = RepairableRBD(
... [("s", "v1"), ("s", "v2"), ("v1", "t"), ("v2", "t")],
... {"v1": valve(), "v2": valve()},
... )
>>> plan = pair.optimal_inspection_intervals(
... allowed=calendar, min_availability=1 - 1e-3
... )
>>> plan.intervals
{'v1': 17520.0, 'v2': 17520.0}
>>> round(1 - plan.availability, 6)
0.000399
The redundant pair meets the target with tests every two years; the
single valve needs them monthly. With a 10% common cause between
the two, and a PFDavg of at most 5e-4, testing them at
different times as well finds a common-cause failure sooner, and
meets the target with yearly tests six months apart, for a third
less:
>>> from repyability import BetaFactor, CCFGroup
>>> common = RepairableRBD(
... [("s", "v1"), ("s", "v2"), ("v1", "t"), ("v2", "t")],
... {"v1": valve(), "v2": valve()},
... ccf_groups=[CCFGroup(["v1", "v2"], BetaFactor(0.1))],
... )
>>> target = 1 - 5e-4
>>> plan = common.optimal_inspection_intervals(
... allowed=calendar, min_availability=target
... )
>>> plan.intervals, round(plan.cost_rate, 4)
({'v1': 4380.0, 'v2': 8760.0}, 0.1712)
>>> staggered = common.optimal_inspection_intervals(
... allowed=calendar,
... min_availability=target,
... offset_shares="stagger",
... )
>>> staggered.intervals, staggered.offsets
({'v1': 8760.0, 'v2': 8760.0}, {'v1': 0.0, 'v2': 4380.0})
>>> round(staggered.cost_rate, 4), round(1 - staggered.availability, 7)
(0.1142, 0.0004925)
with_intervals ¶
with_intervals(intervals, offsets=None)
A copy of this RBD with some maintenance or test intervals changed.
To simulate a plan that optimal_replacement_intervals or
optimal_inspection_intervals chose, with the repair crews it was
chosen without (assume_unlimited_crews), or to compare it
with the schedules this RBD has. The copy is built as this RBD was,
with each named component's schedule given the new interval.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
intervals
|
MaintenancePlan or dict
|
A plan (its |
required |
offsets
|
dict
|
|
None
|
Returns:
| Type | Description |
|---|---|
RepairableRBD
|
The new RBD. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a node has no maintenance or test schedule (or an offset no
test schedule), a block-replaced component is given |
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> pump = {
... "reliability": surv.Weibull.from_params([1000, 2.5]),
... "repairability": surv.Exponential.from_params([0.05]),
... "replace_cost": 5000.0,
... "preventive": {"interval": 1000.0, "cost": 1000.0},
... }
>>> rbd = RepairableRBD([("s", "p"), ("p", "t")], {"p": pump})
>>> plan = rbd.optimal_replacement_intervals()
>>> rbd.with_intervals(plan).expected_cost_rate() == plan.cost_rate
True
initialize_event_queue ¶
initialize_event_queue(t_simulation, working_nodes=None, broken_nodes=None, method='p', sources=None, state=None)
Start one simulation of the system over [0, t_simulation).
Advanced, event-stepping API: availability runs its simulations
through it, and a parent RBD uses it to drive a nested
RepairableRBD, which therefore has the same event API as a
NonRepairable component. Every component starts working and new
(or from its state) except the broken_nodes, which stay down
for the whole window; the working_nodes never fail. The first
failure of each other component is queued (a nested RBD starts its
own simulation and queues its first state change), and events at or
after t_simulation are dropped. Then call next_event
repeatedly to step through the system's state changes; see it for
an example.
Unlike availability, this takes no seed and does not check the
working/broken nodes: the draws come from numpy's global RNG, so
seed that (np.random.seed) for a reproducible run. The state is
kept on the RBD itself (system_state, component_status,
t_simulation, last_change_planned), and availability
uses and then deletes the same state, so do not interleave the two.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t_simulation
|
float
|
The end of the simulated window. |
required |
working_nodes
|
Collection[Hashable]
|
Nodes held working for the whole window, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed for the whole window, by default None. |
None
|
method
|
str
|
Evaluate the system state from the minimal path sets ( |
'p'
|
sources
|
dict
|
Internal: node name -> the object each component's events are
drawn from, which |
None
|
state
|
dict
|
Start the components it names from their states rather than new,
as for |
None
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
If |
mean_unavailability ¶
mean_unavailability(working_nodes=None, broken_nodes=None, method='p')
Returns the system's long-run (steady-state) unavailability.
1 - mean_availability(): the long-run fraction of time the
system is down, exact as
mean_availability
is, and over the same cases. It is worked out in its own right,
from each component's unavailability (MTTR / (MTTF + MTTR),
say, or 1 - exp(-lambda * u) at a time u since a test), as
a sum of products over the structure, so that a small one keeps its
full relative precision: a PFDavg of 1e-12 gets 1e-12,
where one less the availability would keep only about four digits
of it.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working (unavailability 0), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed (unavailability 1), by default None. |
None
|
method
|
str
|
|
'p'
|
Returns:
| Type | Description |
|---|---|
float
|
Long-run unavailability of the system, in |
Raises:
| Type | Description |
|---|---|
ValueError
|
As for |
NotImplementedError
|
As for |
Examples:
Two components in parallel, each down a millionth of the time, are down together a millionth of a millionth of it:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> unit = {
... "reliability": surv.Exponential.from_params([1e-6]),
... "repairability": surv.ExactEventTime.from_params([1.0]),
... }
>>> rbd = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {"a": unit, "b": unit},
... )
>>> f"{rbd.mean_unavailability():.6e}"
'9.999980e-13'
analysis_routes ¶
analysis_routes()
How each analysis of this RBD is computed, found without running it: exactly, numerically, by simulation, or not at all.
The route follows from the method, then from the components: the
long-run values are exact from each component's long-run
availability, which is a closed form, or numerical for a component
under preventive maintenance, or refused for a case they do not
cover (hidden failures at a rate that is not constant, say); the
availability over time solves each component's renewal equation on
a grid; availability, cost and compare always simulate.
A component whose life has no exact or numerical reliability (a
StandbyModel or LoadSharingModel that only simulations take)
makes the values it enters refuse, as it does. A refusal is found by
the check the method itself runs, and its reason is the message it
would raise; limits that only computing shows (a grid grown too
large) are not foreseen.
A method that takes nodes, targets or intervals is described for
its defaults: optimal_replacement_intervals choosing for every
component under age replacement, say.
For a simulation, engine is the engine engine="auto" runs a
long simulation on: when the compiled engine simulates the system,
"numba" if numba is installed (or the name of an engine another
package adds, see repyability.rbd.engines), else "python",
with the reason in engine_reason.
Returns:
| Type | Description |
|---|---|
dict[str, AnalysisRoute]
|
For each public analysis, by method name: its |
Examples:
A pump found failed only by monthly tests that take it off line for about an hour: its values are numerical, its cycle from one test that finds a failure to the next followed test by test. Tests that could last as long as the month (of 1000 hours on average, say) are refused, and the simulation is the way:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def pump(test_hours):
... test = surv.Exponential.from_params([1.0 / test_hours])
... return RepairableRBD(
... [("s", "pump"), ("pump", "t")],
... {
... "pump": {
... "reliability": surv.Weibull.from_params([500, 2]),
... "repairability": "instant",
... "inspection": {"interval": 720, "duration": test},
... }
... },
... )
>>> pump(1.0).analysis_routes()["mean_availability"].route
'numerical'
>>> routes = pump(1000.0).analysis_routes()
>>> routes["mean_availability"].route
'refused'
>>> routes["mean_availability"].nodes
('pump',)
>>> routes["availability"].route
'simulated'
mean_availability ¶
mean_availability(working_nodes=None, broken_nodes=None, method='p')
Returns the system's long-run (steady-state) availability.
Exact, with no simulation. Each component's long-run availability
is MTTF / (MTTF + MTTR) (1 if it is repaired instantly; a nested
RepairableRBD's is its own mean_availability; under age
replacement, see node_availability). Since the components fail
and are repaired independently, the system's is
the RBD's structure function evaluated exactly at those
availabilities. It is the long-run fraction of time the system is
up; for availability over time, from a start with everything
working, use availability.
A component with hidden failures is up with probability
exp(-lambda * u) at a time u since its last inspection, so
components inspected at the same times are down together more
often than independent ones would be. The system's availability is
then averaged over time, over one period of the inspection
schedules (the least common multiple of their intervals): for two
such components in parallel, inspected together every tau, the
unavailability is (1 / tau) * integral_0^tau (1 -
exp(-lambda * t)) ** 2 dt, about (lambda * tau) ** 2 / 3: the
average probability of failure on demand (PFDavg) of a 1oo2 safety
function. This is exact for a constant failure rate, and numerical
(summed over the test intervals) for any other life, with instant
tests and instant repair (see node_availability).
With fewer repair_crews than components, a failed component can
wait for a crew, and the components no longer fail and recover
independently. When the components the crews work on have
exponential lives and exponential (or instant) repairs, with no
scheduled maintenance or inspection, the system is a Markov chain:
its state is which components are under repair and which are
waiting, in the order the crews will take them. Its long-run
distribution is solved exactly (for up to 15,000 states), and the
availability is the system's, state by state, averaged over it. A
node held working never fails and one held broken is never
repaired, so neither needs a crew: the chain is of the others. A
nested RepairableRBD has crews of its own, and enters through
its own long-run availability.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working (availability 1), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed (availability 0), by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
Returns:
| Type | Description |
|---|---|
float
|
Long-run availability of the system, in |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
If a component's life has no exact or numerical mean (a
|
Examples:
A single component with mean time to failure 10 and mean time to
repair 1 has long-run availability 10 / (10 + 1):
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
... },
... )
>>> round(rbd.mean_availability(), 4)
0.9091
point_availability ¶
point_availability(x, working_nodes=None, broken_nodes=None, method='p', state=None)
The probability that the system is up at each time x, with
every component new at 0: exact, with no simulation.
Each component's point availability A(t) follows from the
distributions of its up and down times by the renewal equation,
solved numerically (see repyability/rbd/_point_availability.py):
its error is about 4e-7 (up to 4e-6 soon after the start). Since
the components fail and are repaired independently, the system's is
the structure function evaluated exactly at theirs, at each time.
It starts at 1 (less any components dead on arrival) and settles at
mean_availability
(or, with components replaced or inspected on a calendar, repeats
with the calendar about it);
availability estimates
the same curve by simulation. At a time something happens -- a
block replacement, say -- it is the availability just after.
Components under age or block replacement are covered, as in
mean_availability, and nested RBDs through their own point
availability. A component with hidden failures is up with the
probability that it has not failed since its last test, u
before (exp(-lambda * u) for a constant failure rate), or
after it with its tests and repairs taking time or its tests
missing some failures: as in mean_availability.
With fewer repair_crews than components, a component can wait
for a crew and the components are not independent: the system's
point availability is then that of the crews' Markov chain (see
mean_availability), followed over time from its state at 0 by
uniformization, to about 1e-13; a nested RBD, which has crews of
its own, enters through its own point availability. A standby
group's comes from its units' chain in the same way, and a
common-cause group's members' joint states from theirs (see
ccf_groups), followed from every member up through their tests
or by uniformization, with the system summed over them.
Each component's curve is computed on a grid of 1,000 steps over its
typical up time. Near a time at which its units start or stop on a
schedule -- at 0, at its scheduled replacements, and at the failures
of a lifetime known exactly -- what happens faster than a step, such
as a short repair, is smoothed over the step: a point value within
a step of such a time can be off by up to about the probability
that the component is under repair then. mission_availability
is not affected.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array - like
|
Times, from 0 (every component new). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work (availability 1), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed (availability 0), by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
Returns:
| Type | Description |
|---|---|
float or ndarray
|
The system's point availability at each time, in |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a time is negative or not finite, or a working/broken node is
invalid (see |
NotImplementedError
|
If a component has hidden failures its numerical values do not
cover, a model that block replacement's exact values do not
cover (see |
Examples:
One component with failure rate 0.1 and repair rate 1 is up at time
t with probability 1 / 1.1 + (0.1 / 1.1) * exp(-1.1 t):
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
... },
... )
>>> rbd.point_availability([0.0, 1.0, 100.0]).round(4).tolist()
[1.0, 0.9394, 0.9091]
point_unavailability ¶
point_unavailability(x, working_nodes=None, broken_nodes=None, method='p', state=None)
The probability that the system is down at each time x, with
every component new at 0: one less point_availability, worked
out in its own right so that a small one keeps its precision,
as mean_unavailability is in the long run.
Each component's probability of being down enters the structure
function's sum of products for the system's failing, so that a
redundant system's 1e-16 is not lost below one less a number
next to 1, as 1 - point_availability(t) loses it; a
common-cause group's members' joint states enter theirs in the same
way. A component new at 0 with an exponential life and an
exponential (or instant) repair, and no maintenance or tests, is
down with lambda / (lambda + mu) (1 - exp(-(lambda + mu) t)),
in closed form; any other, with one less its point availability
(see
point_availability),
whose grid smooths over a step what changes faster than it. With
components waiting for repair crews, the crews' chain gives it as
one less the availability (to about 1e-13).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array - like
|
Times, from 0 (every component new). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work (never down), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed (always down), by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new,
as for |
None
|
Returns:
| Type | Description |
|---|---|
float or ndarray
|
The system's point unavailability at each time, in |
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
Examples:
Two valves in parallel, each failing at a rate of 1e-9 an hour
and repaired in a mean of eight hours: each is down 8e-9 of
the time, and the pair, both at once, 6.4e-17, which one less
the availability rounds to 0 or 1.1e-16:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> valve = {
... "reliability": surv.Exponential.from_params([1e-9]),
... "repairability": surv.Exponential.from_params([1 / 8]),
... }
>>> pair = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {"a": valve, "b": valve},
... )
>>> f"{pair.point_unavailability(1000.0):.3e}"
'6.400e-17'
mission_availability ¶
mission_availability(t, working_nodes=None, broken_nodes=None, method='p', state=None)
The expected fraction of [0, t] the system is up, with every
component new at 0: exact, with no simulation.
It is the mean of
point_availability
over the window, integrated by 4-point Gauss-Legendre quadrature on
pieces that end where the components' curves bend (a scheduled
replacement, a test, a down time from a known instant) and are at
most a few steps of the finest grid the curves are worked out on,
each halved until the quadrature on it agrees with that on its
halves, to about 1e-8: it takes little longer than the point
availability. That runs up to the time when the curves have all
settled at their long-run values (or into repeating with their
inspections and block replacements); past it, the integral is
extended exactly, so a mission of decades costs no more than one of
a few years. It is what
availability estimates
by simulation as the mean of each simulation's uptime divided by
t_simulation. As t grows it approaches mean_availability,
from which it differs by about b / t, b a constant of the
components' up and down times: positive unless their lives vary
more than exponential ones do, and largest for components that wear
out (Weibull shape above 1, say), which from new fail less early
on.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The windows' lengths (one mission average for each). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work (availability 1), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed (availability 0), by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
Returns:
| Type | Description |
|---|---|
float or ndarray
|
The mission availability for each window, in |
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
Examples:
One component with failure rate 0.1 and repair rate 1, over 10 time
units: 1 / 1.1 + 0.1 / 1.1 ** 2 * (1 - exp(-11)) / 10.
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
... },
... )
>>> round(rbd.mission_availability(10.0), 4)
0.9174
mission_unavailability ¶
mission_unavailability(t, working_nodes=None, broken_nodes=None, method='p', state=None)
The expected fraction of [0, t] the system is down, with
every component new at 0: one less mission_availability, worked
out in its own right so that a small one keeps its precision:
the mean of
point_unavailability
over the window, integrated as mission_availability integrates
the availability.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The missions' lengths, from 0 (every component new). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work (never down), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed (always down), by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new,
as for |
None
|
Returns:
| Type | Description |
|---|---|
float or ndarray
|
The expected fraction of each mission the system is down, in
|
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
Examples:
The pair of valves of point_unavailability, down together
6.4e-17 of the time once they settle, from new over a year:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> valve = {
... "reliability": surv.Exponential.from_params([1e-9]),
... "repairability": surv.Exponential.from_params([1 / 8]),
... }
>>> pair = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {"a": valve, "b": valve},
... )
>>> f"{pair.mission_unavailability(8760.0):.3e}"
'6.391e-17'
expected_failures ¶
expected_failures(t, working_nodes=None, broken_nodes=None, method='p', state=None)
The expected number of system failures in [0, t), with every
component new at 0: exact, with no simulation.
A component's failure takes the system down if the component is
critical then: if the system is up with it and down without it.
The components fail and are repaired independently, so at a time
s that is so with probability I_B^i(s), the component's
Birnbaum importance at the components' point availabilities then,
and the system's expected failures are
integral from 0 to t of sum_i I_B^i(s) dM_i(s)
M_i(s) the component's expected failures by s: the
time-dependent form of the Birnbaum/Vesely formula, whose long-run
rate is system_failure_frequency. Each M_i follows from the
component's renewal equation, solved on the same grid as its point
availability (see point_availability), to about 1e-7; the
integral is summed over the pieces mission_availability
integrates over, on each with the importance taken as the cubic
through its values at the quadrature points, and extended exactly
past the time the curves have settled. Failures at exact
times (units dead on arrival, an exact lifetime) are taken together
when several fall at once. It is what availability estimates as
system_failures / n_simulations: every system failure counts,
including the zero-length outages an instantly repaired component
causes; planned outages do not (see expected_events).
With fewer repair_crews than components, the components are not
independent: the rate of system failures is then that of the crews'
Markov chain in each of its states (a component up failing at its
rate where it is critical), integrated over the chain's states over
time by uniformization, to about 1e-13 (see point_availability).
The crews' chain takes no nested RBD's events in, as yet.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The windows' lengths (one count for each). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work, and never fail, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed, by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
Returns:
| Type | Description |
|---|---|
float or ndarray
|
The expected number of system failures in each window, in
|
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
Examples:
One component with failure rate 0.1 and repair rate 1 fails at the
rate 0.1 while it is up, 0.1 * (10 / 1.1 + 0.1 / 1.1 ** 2 * (1 -
exp(-11))) times in 10 time units:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
... },
... )
>>> round(rbd.expected_failures(10.0), 4)
0.9174
expected_events ¶
expected_events(t, working_nodes=None, broken_nodes=None, method='p', state=None)
What the system and each component are expected to do in
[0, t), with every component new at 0: exact, with no
simulation.
The system's failures are those of expected_failures, and its
planned outages are counted the same way, from its components'
preventive maintenance that takes time; its down time is t * (1
- mission_availability(t)). Each component's failures, corrective
and preventive actions and tests follow from its renewal equation,
solved on the grid of its point availability (see
point_availability), to about 1e-7, and its down time is the
integral of its unavailability. The maintenance due at exact times
(at a replacement age, at block times) is counted exactly, and
maintenance of several components due at the same time takes the
system down at most once. It is the mean of what availability
counts in each simulation of the window, and events at t itself
fall outside the window, as there.
A component's corrective actions, at each of which its
repair_cost and replace_cost are charged, are its failures,
except for hidden failures: those found by a test in the window. A
nested RBD's failures are its system failures, and its own
maintenance is not this RBD's to count. A forced node does nothing:
held working it is never down, held broken it is down throughout.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The windows' lengths. |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed, by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
Returns:
| Type | Description |
|---|---|
ExpectedEvents
|
The system's expected failures, planned outages and down time,
and each component's expected failures, corrective and
preventive actions, tests and down time: floats for one window,
arrays in |
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
Examples:
A pump that wears out, replaced at 600 hours in about 7 (a planned outage), over a year:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "p"), ("p", "t")],
... {
... "p": {
... "reliability": surv.Weibull.from_params([1000, 2.5]),
... "repairability": surv.LogNormal.from_params([3, 0.5]),
... "preventive": {
... "interval": 600.0,
... "duration": surv.Weibull.from_params([8, 3]),
... },
... }
... },
... )
>>> year = rbd.expected_events(8760.0)
>>> round(year.system_failures, 3), round(year.node_preventive["p"], 3)
(3.706, 11.275)
>>> round(year.system_planned_outages, 3) # each replacement
11.275
expected_cost ¶
expected_cost(t, working_nodes=None, broken_nodes=None, method='p', state=None, *, discount_rate=0.0)
The expected cost of running the system for t, from new:
exact, with no simulation.
Each category is its events' expected number over [0, t) (see
expected_events) times its mean cost, as cost charges them:
repair_cost and replace_cost at each corrective action, the
preventive cost at each preventive replacement, the inspection cost
at each test, downtime_cost over each component's expected down
time, downtime_cost_rate over the system's, and a maintenance
group's set-up cost once per stop (each failure or replacement of a
member, those at one instant one stop). It is the mean that
cost estimates by simulation, with its distribution;
total_cost is the long-run approximation, exact only over a
window long compared with the components' cycles. As t grows
this approaches expected_cost_rate() * t plus a constant, of the
early years from new.
Every cost is optional, and a cost given as a distribution enters
through its mean. Only this RBD's own costs count, not those inside
a nested RepairableRBD. With nothing priced, not even an
acquisition cost (see has_costs), there is no cost model to
evaluate, and this refuses; with an acquisition cost alone every
category is 0, exactly.
With a discount_rate r each value is its present value,
the costs discounted continuously from when they fall:
integral from 0 to t of exp(-r s) dC(s), C the expected
cost from new, worked out by parts, exp(-r t) C(t) + r *
integral from 0 to t of exp(-r s) C(s) ds, the integral by
Gauss-Kronrod quadrature on pieces halved until the total is
within 1e-8 of itself. It is the exact form of total_cost's
discounted value, which spends the long-run rate from the start:
the early years, which count most, are those that differ from it.
The acquisition cost, at the start, is not discounted.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The windows' lengths. |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that always work: they incur no corrective, preventive, inspection or downtime cost, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always failed: they incur no corrective cost, but their own downtime cost throughout, by default None. |
None
|
method
|
str
|
Evaluate the structure function from the minimal path sets
( |
'p'
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
discount_rate
|
float
|
The continuous discount rate per unit time of the component
models, as for |
0.0
|
Returns:
| Type | Description |
|---|---|
ExpectedCost
|
The expected cost of each window, by category and component,
with the acquisition cost beside it: present values with a
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If nothing is priced, not even an acquisition cost (see
|
NotImplementedError
|
As for |
Examples:
A pump bought for 20,000, failing on average every 1000 hours, repaired in 10 at 500 per repair, with lost production at 100 per hour, owned for a year:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "p"), ("p", "t")],
... {
... "p": {
... "reliability": surv.Exponential.from_params([1e-3]),
... "repairability": surv.Exponential.from_params([0.1]),
... "repair_cost": 500.0,
... "acquisition_cost": 20000.0,
... }
... },
... downtime_cost_rate=100.0,
... )
>>> year = rbd.expected_cost(8760.0)
>>> round(year.mean), round(year.total)
(13000, 33000)
>>> round(rbd.total_cost(8760.0)) # at the long-run rate
33010
Owned for ten years at 7% a year, from new and at the long-run rate:
>>> import math
>>> r = math.log(1.07) / 8760 # per hour
>>> round(rbd.expected_cost(87600.0, discount_rate=r).total)
114528
>>> round(rbd.total_cost(87600.0, discount_rate=r))
114538
capacity_distribution ¶
capacity_distribution(working_nodes=None, broken_nodes=None)
The exact long-run distribution of the system's capacity.
Each node carries its capacity (given when the RBD was built)
while it is up, at one level or several, and nothing while it is
down; the system's capacity is the most that can flow through the
nodes that are up from the input to the output. A node given no
capacity limits nothing, unless its model has capacities of its own:
a degrading component spends its up time in its stages in
proportion to their mean times (by the renewal-reward theorem), and
a nested RBD with capacities brings its own long-run
distribution. The probability of each level is the long-run fraction of
time the system spends at it; the capacity is positive exactly when
the system is up, so the fraction of time it is positive is
mean_availability.
With three pumps of half the demand each, one down costs nothing
and two cost half the output.
Exact, with no simulation: worked out from each node's long-run
availability as mean_availability is (see
RBD.system_capacity), and
averaged over the schedules of components inspected or replaced on
a calendar, which are down together more often than independent
ones would be. With limited repair_crews, it is averaged over
the states of the crews' Markov chain instead (see
mean_availability).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always being up, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being down, by default None. |
None
|
Returns:
| Type | Description |
|---|---|
CapacityDistribution
|
The capacities the system can have, and the long-run fraction of
time at each. Its |
Raises:
| Type | Description |
|---|---|
ValueError
|
If no node has a capacity, or as for |
NotImplementedError
|
As for |
Examples:
Three pumps of 50 each, each up 10 / 11 of the time, against a demand of 100:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> pump = {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
>>> plant = RepairableRBD(
... [("s", "a"), ("s", "b"), ("s", "c"),
... ("a", "t"), ("b", "t"), ("c", "t")],
... {"a": pump, "b": pump, "c": pump},
... capacity={"a": 50, "b": 50, "c": 50},
... )
>>> capacity = plant.capacity_distribution()
>>> capacity.levels.tolist()
[0.0, 50.0, 100.0, 150.0]
>>> round(capacity.meets(100), 4) # at least two up
0.9767
>>> round(capacity.delivered_fraction(100), 4)
0.988
>>> round(capacity.mean, 2) # 150 * 10 / 11
136.36
point_capacity ¶
point_capacity(x, working_nodes=None, broken_nodes=None, state=None)
The distribution of the system's capacity at each time x,
with every component new at 0: exact, with no simulation.
It is capacity_distribution at a time rather than in the long
run: the components fail and are repaired independently, so each
is up at t with its point availability (see
point_availability), and the system's capacity distribution is
the same exact computation at those (see
RBD.system_capacity). A node
with levels while it works is at each with the level's share of its
availability. A degrading component (a DegradingNode as its
reliability) is in each stage with the probability that its first
unit is, or a unit put into service later is, by its renewals,
solved on the grid of its availability; a nested RBD with
capacities brings its own distribution over time. The probability
of a capacity above 0 is point_availability(x), and the
distribution settles at capacity_distribution();
availability(demand=...) estimates its mean by simulation.
With fewer repair_crews than components the components are not
independent: the probability of each capacity is then that of the
states of the crews' Markov chain at it, followed over time by
uniformization (see point_availability). The chain takes no
nested RBD's capacities in, as yet.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array - like
|
Times, from 0 (every component new). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that are always up, by default None. One that takes its capacity from its model (a degrading component, a nested RBD) is at the levels it is up at, in proportion. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always down, by default None. |
None
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
Returns:
| Type | Description |
|---|---|
CapacityDistribution
|
The capacities the system can have and their probabilities:
one per level for a scalar |
Raises:
| Type | Description |
|---|---|
ValueError
|
If no node has a capacity, or as for |
NotImplementedError
|
If a degrading component is maintained or inspected on a
schedule, or as for |
Examples:
Three pumps of 50 each, new at 0, with failure rate 0.1 and repair
rate 1, each up at t with probability 1 / 1.1 + 0.1 / 1.1 *
exp(-1.1 t), against a demand of 100:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> pump = {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
>>> plant = RepairableRBD(
... [("s", "a"), ("s", "b"), ("s", "c"),
... ("a", "t"), ("b", "t"), ("c", "t")],
... {"a": pump, "b": pump, "c": pump},
... capacity={"a": 50, "b": 50, "c": 50},
... )
>>> capacity = plant.point_capacity([0.0, 1.0, 100.0])
>>> capacity.levels.tolist()
[0.0, 50.0, 100.0, 150.0]
>>> capacity.meets(100).round(4).tolist() # at least two up
[1.0, 0.9894, 0.9767]
>>> round(plant.capacity_distribution().meets(100), 4) # long run
0.9767
mission_capacity ¶
mission_capacity(t, working_nodes=None, broken_nodes=None, state=None)
The distribution of the system's capacity over [0, t], with
every component new at 0: the expected fraction of the window it
spends at each level. Exact, with no simulation.
It is the mean of point_capacity over the window, integrated as
mission_availability integrates the availability: by
Gauss-Legendre quadrature on pieces cut where the components' curves
bend, and extended exactly past the time they have settled.
Its delivered_fraction(demand) is the expected fraction of the
demand delivered over the window, the production availability of a
contract period from new, which availability(demand=...)
estimates by simulation as delivered_fraction; meets(demand)
is the expected fraction of the window the capacity meets the
demand, and mean() the average capacity. As t grows it
approaches capacity_distribution().
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The windows' lengths (one distribution for each). |
required |
working_nodes
|
Collection[Hashable]
|
Nodes that are always up, by default None (see
|
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes that are always down, by default None. |
None
|
state
|
dict or str
|
Start from the components' current states rather than new:
|
None
|
Returns:
| Type | Description |
|---|---|
CapacityDistribution
|
The capacities and the expected fraction of each window spent
at each: one per level for a scalar |
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
Examples:
The three pumps above over 10 time units, against a demand of 100: more than in the long run, as every pump starts up.
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> pump = {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
>>> plant = RepairableRBD(
... [("s", "a"), ("s", "b"), ("s", "c"),
... ("a", "t"), ("b", "t"), ("c", "t")],
... {"a": pump, "b": pump, "c": pump},
... capacity={"a": 50, "b": 50, "c": 50},
... )
>>> round(plant.mission_capacity(10.0).delivered_fraction(100), 4)
0.9896
>>> round(plant.capacity_distribution().delivered_fraction(100), 4)
0.988
next_event ¶
next_event(method='p', sources=None)
Advance the current simulation to the system's next state change.
Advanced, event-stepping API, after initialize_event_queue. It
takes the queued component events in time order, updating each
component's state and queueing its next event (if before
t_simulation), until the system changes state, and returns
that change. It gives a RepairableRBD the same event API as a
NonRepairable, so it can be a node of another RBD. Unlike
NonRepairable.next_event, which returns the time to the
component's next event, this returns the time of the change, on
the simulation's clock.
When no further change happens before t_simulation, it returns
(t_simulation, current state) and discards the queue, so the
next call raises a ValueError until initialize_event_queue is
called again. After each call last_change_planned says whether
the change was the start of a planned outage (preventive
maintenance that takes the system down).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
method
|
str
|
Evaluate the system state from the minimal path sets ( |
'p'
|
sources
|
dict
|
Internal: the same |
None
|
Returns:
| Type | Description |
|---|---|
tuple[float, bool]
|
The time of the change, and the system's new state (True for working, False for failed). |
Raises:
| Type | Description |
|---|---|
ValueError
|
If the event queue has not been initialised, or was used up by
an earlier call; or if |
Examples:
A component that fails exactly 10 hours after each repair, and takes exactly 2 hours to repair:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.ExactEventTime.from_params(10),
... "repairability": surv.ExactEventTime.from_params(2),
... }
... },
... )
>>> rbd.initialize_event_queue(30.0)
>>> changes = [rbd.next_event()]
>>> while changes[-1][0] < 30.0:
... changes.append(rbd.next_event())
>>> changes[:4]
[(10.0, False), (12.0, True), (22.0, False), (24.0, True)]
The next failure, at 34, is outside the window, so the last call returns the end of the window and the state the system is left in:
>>> changes[4]
(30.0, True)
availability ¶
availability(t_simulation, working_nodes=None, broken_nodes=None, method='p', mc_samples=None, verbose=False, seed=None, *, tolerance=None, confidence=0.95, max_samples=None, antithetic=False, n_jobs=None, demand=None, engine='auto', state=None, curve_points=None, shard_map=None, shard_size=None, control_variate=None, conditional=None)
Simulate the system's availability over [0, t_simulation].
Runs mc_samples independent Monte-Carlo (discrete-event)
simulations of the system from time 0, each starting with every
component working and new (except broken_nodes), or from the
components' current state. Each component alternates failure and
repair independently of the others and of the system state: it
fails after a time drawn from its reliability model and is
restored, as good as new, after a time drawn from its repairability
model. A component under preventive maintenance is also replaced on
its schedule (see the class docstring), down for the maintenance
time if it takes any: a planned outage, which counts in every up
and down time but is not a failure. A nested RepairableRBD runs
its own simulation on the same clock, and is down while its system
is down. The system state is re-evaluated at every component event.
The result holds the estimated availability over time, the system's
and every node's total up and down time, the number of system
failures and restorations, and criticality measures that attribute
the system's up and down time, failures and restorations to the
nodes (see Criticalities). When
something is priced (see has_costs), the same simulations also
accumulate costs, and the result's cost is a
CostResult; otherwise it is None. For
exact long-run values, use mean_availability and the other
steady-state methods.
When nodes have capacities (see the capacity argument of the
class), the same simulations also follow what the system can
deliver: its capacity after every component event (a failure that
leaves the system up can still take some of its capacity), the time
it spends at each capacity, and the fraction of demand it
delivers (see
AvailabilityResult). The
simulation draws only whether each component is up, so a node
working at several levels counts at each in proportion to its
probability. Nodes taking their capacity from their models (a
DegradingNode's stages, or a nested RBD's capacities) are not
followed: give them a capacity, or use capacity_distribution for
the long run.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t_simulation
|
float
|
Length of the window each simulation covers, in the time unit of the component models. Events at or after it are not simulated. |
required |
working_nodes
|
Collection[Hashable]
|
Nodes held working for the whole window: they never fail. By default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed for the whole window: they are down from time 0 and never repaired. By default None. |
None
|
method
|
str
|
Evaluate the system state from the minimal path sets ( |
'p'
|
mc_samples
|
int
|
Number of simulations, by default 10_000. |
None
|
verbose
|
bool
|
If True, displays a progress bar of the simulations, by default False. |
False
|
seed
|
int
|
Seed for a reproducible run, by default None. Each random
quantity -- a component's times to failure, its repair times,
its maintenance or test times, a cost given as a distribution
-- is drawn from a stream of its own, named by the component's
place and the quantity, and seeded from |
None
|
tolerance
|
float
|
Simulate until the mean availability over the window (the
fraction of it the system is up) is known to within
|
None
|
confidence
|
float
|
The confidence level |
0.95
|
max_samples
|
int
|
The most simulations a run to |
None
|
antithetic
|
bool
|
Run the simulations in antithetic pairs, by default False: in
the second simulation of a pair, each stream (see |
False
|
n_jobs
|
int
|
Run the simulations in parallel on |
None
|
engine
|
str
|
What runs the simulations: |
'auto'
|
demand
|
float
|
The demand the delivered fraction is measured against, in the capacities' units, when nodes have capacities. By default the system's capacity with every component up (at its highest level): its design capacity. If that is unlimited, no delivered fraction is worked out unless a demand is given. |
None
|
state
|
dict
|
Start each simulation from the components' current states
rather than new: |
None
|
curve_points
|
int
|
Keep the availability over time on a grid of |
None
|
shard_map
|
callable
|
Run the simulations as shards (see |
None
|
shard_size
|
int
|
The simulations to a shard, with |
None
|
control_variate
|
bool
|
Control the estimate of the mean availability by the system's
exact twin: the same diagram, its components failing and
repaired independently, whose mean availability over the
window is exact, simulated alongside with common random
numbers. By default None: where the exact methods work out the
system's expected values over the window, the result's means
( |
None
|
conditional
|
bool
|
Take the expected values given the histories of the dependent
modules: the nodes whose values over time the exact
methods do not work out ( |
None
|
Returns:
| Type | Description |
|---|---|
AvailabilityResult
|
The availability over time ( |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output
node, or is in both sets, if |
NotImplementedError
|
With |
ImportError
|
With |
Examples:
Two components in series, each with MTTF 10 and MTTR 1:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> unit = {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... }
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")], {"a": unit, "b": unit}
... )
>>> result = rbd.availability(t_simulation=50, mc_samples=200, seed=0)
>>> float(result.timeline[0]), float(result.availability[0])
(0.0, 1.0)
>>> float(result.timeline[-1])
50.0
The fraction of the window the system was up is close to the
long-run availability, (10 / 11) ** 2 (over a short window it is
a little above it on average, as every simulation starts up):
>>> window = result.n_simulations * result.time_simulated_to
>>> round(float(result.system_uptime) / window, 4)
0.832
>>> round(rbd.mean_availability(), 4)
0.8264
In series the system is up only while every node is up:
>>> oci = result.criticalities.operational_criticality_index
>>> {node: round(float(v), 4) for node, v in oci.up.items()}
{'a': 1.0, 'b': 1.0}
simulate_timelines ¶
simulate_timelines(t_simulation, mc_samples=None, seed=None, working_nodes=None, broken_nodes=None, antithetic=False, engine='auto', n_jobs=None, start=0, state=None)
Simulate the system over [0, t_simulation] and keep each
simulation's up/down histories as timelines: every component's and
the system's, with the component that caused each of the system's
changes. A common-cause group's shared cause fails the members it
strikes one after another at one instant, each a failure of the
member's own, so a system failure it makes is the member's whose
fall took the system down (the last the structure could spare),
not the group's: failures_by_cause counts it under that
member.
These are the simulations
availability runs with
the same seed (and options), kept whole rather than added up:
so any measure of a history can be read off them (the time to the
first system failure, the longest outage, the outages per year),
and they merge with other timelines (see
repyability.timelines). Their histories are the event loop's,
whichever engine makes them: each simulation's system uptime is
availability's to the last bit.
The event loop records them as it runs, compiled with numba where
it can (pip install "repyability[fast]"), as availability
runs. On the Python engine, independent components (plain units
with streamed models, and nested RBDs of them) have their histories
drawn straight from their streams instead, a batch of simulations
at once, and the system's merged from theirs (see
RBD.system_timeline): a
simulation in which components change at the same instant is run
in the event loop, which orders them as it does.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t_simulation
|
float
|
The window's end. |
required |
mc_samples
|
int
|
The number of simulations. By default 1,000. |
None
|
seed
|
int
|
Seeds the simulations, as for |
None
|
working_nodes
|
collection
|
Nodes held working throughout. |
None
|
broken_nodes
|
collection
|
Nodes held broken throughout. |
None
|
antithetic
|
bool
|
Simulate in antithetic pairs ( |
False
|
engine
|
str
|
|
'auto'
|
n_jobs
|
int
|
Run on this many processes (the Python event loop) or threads
(numba, and the streams), |
None
|
start
|
int
|
The run's first simulation to make: simulations |
0
|
state
|
dict
|
Start from the components' current states rather than new, as
for |
None
|
Returns:
| Type | Description |
|---|---|
TimelineSimulation
|
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
If a common-cause group's members are not exponential, are
held, or start from a state (see |
ImportError
|
With |
Examples:
Two pumps in parallel: how long until the plant first fails, and which pump's failure did it:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> unit = {
... "reliability": surv.Weibull.from_params([100, 1.5]),
... "repairability": surv.Exponential.from_params([0.5]),
... }
>>> rbd = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {"a": unit, "b": unit},
... )
>>> runs = rbd.simulate_timelines(1000.0, mc_samples=2000, seed=1)
>>> bool(runs.system.first_failure.min() > 0)
True
>>> sorted(runs.system.failures_by_cause())
['a', 'b']
The same simulations availability runs:
>>> result = rbd.availability(1000.0, mc_samples=2000, seed=1)
>>> bool(np.array_equal(runs.system.uptime, result.uptimes))
True
simulate_chunk ¶
simulate_chunk(t_simulation, start, stop, *, seed, working_nodes=None, broken_nodes=None, method='p', antithetic=False, demand=None, engine='auto', n_jobs=None, verbose=False, state=None, curve_points=None, control_variate=None, conditional=None)
Run simulations start to stop - 1 of the run
availability(t_simulation, mc_samples=N, seed=seed, ...) makes,
for any N of at least stop, and return them as a chunk to
merge with the run's others.
Each simulation draws from streams of its own, seeded from seed
and its position in the run (see
Random streams), so it
comes out the same wherever and whenever it runs, by either engine,
in any company. A run can so be split across processes, machines or
preemptible workers: each runs its chunk, saves it
(SimulationChunk.to_json) and sends it back, and
availability_from_chunks merges the chunks into the run's
AvailabilityResult. Chunks of
simulations 0 to N - 1 give the result of
availability(..., mc_samples=N), to the last bit: the same
simulations, so the same per-simulation values and timeline, and
the same totals, which are kept exactly however the run is cut.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t_simulation
|
float
|
The window each simulation covers, from 0. |
required |
start
|
int
|
The first simulation to run, by its position in the run (from
0; even with |
required |
stop
|
int
|
The simulation after the last to run (above |
required |
seed
|
int
|
The run's seed, required: chunks of a run share it. |
required |
working_nodes
|
Collection[Hashable]
|
Nodes held working, as for |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held broken, as for |
None
|
method
|
str
|
|
'p'
|
antithetic
|
bool
|
Simulate in antithetic pairs, as for |
False
|
demand
|
float
|
A demand on the system's capacity, as for |
None
|
engine
|
str
|
|
'auto'
|
n_jobs
|
int
|
Run the chunk on several CPUs, as for |
None
|
verbose
|
bool
|
Show a progress bar, by default False. |
False
|
state
|
dict
|
The components' states at the start, as for |
None
|
curve_points
|
int
|
Count the curve on a grid, as for |
None
|
control_variate
|
bool
|
None (the default) or False, as for |
None
|
conditional
|
bool
|
None (the default) or False, as for |
None
|
Returns:
| Type | Description |
|---|---|
SimulationChunk
|
The simulations' totals, with the run's settings. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
As for |
See Also
shards : The run's simulations as plain data, to run anywhere.
Examples:
Two chunks of a run, merged, are the run:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "t")],
... {"a": {"reliability": E([0.1]), "repairability": E([1.0])}},
... )
>>> first = rbd.simulate_chunk(100.0, 0, 300, seed=1)
>>> rest = rbd.simulate_chunk(100.0, 300, 500, seed=1)
>>> merged = rbd.availability_from_chunks([first, rest])
>>> whole = rbd.availability(100.0, mc_samples=500, seed=1)
>>> bool((merged.uptimes == whole.uptimes).all())
True
shards ¶
shards(t_simulation, mc_samples=None, *, seed=None, size=None, working_nodes=None, broken_nodes=None, method='p', antithetic=False, demand=None, engine='auto', state=None, curve_points=None)
The run availability(t_simulation, mc_samples, seed=seed,
...) makes, as shards: ranges of its simulations as plain data,
to run anywhere (see repyability.rbd.shards).
Each shard is JSON: this system as to_dict gives it, the run's
entropy (drawn once, here, from the seed) and settings, and its
range. repyability.run_shard runs one in any process, on any
machine and by any executor (concurrent.futures, Ray, Dask, a
batch system with ``python -m repyability.rbd.shards < shard.json
partial.npz
), and gives back its partial, the simulations' totals as.npzbytes;availability_from_chunks(partials, mc_samples)puts them together, in any order, into the run's result, to the last bit.availability(..., shard_map=...)does it all through amap`` of your choosing.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t_simulation
|
float
|
The window each simulation covers, from 0. |
required |
mc_samples
|
int
|
The run's simulations, by default 10_000. |
None
|
seed
|
int
|
The run's seed, by default None: drawn from numpy's global RNG,
as |
None
|
size
|
int
|
Simulations to a shard, rounded up to a whole number of the run's widest block of draws (see Random streams), so that no two shards draw the same block; by default as many as make 1024 or more. A shard should run for some seconds, to repay a worker's start. |
None
|
working_nodes
|
Collection[Hashable]
|
Nodes held working, as for |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed, as for |
None
|
method
|
str
|
|
'p'
|
antithetic
|
bool
|
Run the simulations in antithetic pairs, as for
|
False
|
demand
|
float
|
The demand the delivered fraction is measured against, as for
|
None
|
state
|
dict
|
The components' states at the start, as for |
None
|
engine
|
str
|
What runs each shard: as for |
'auto'
|
curve_points
|
int
|
Count the curve on a grid, as for |
None
|
Returns:
| Type | Description |
|---|---|
list of bytes
|
The shards, in order of their simulations, as JSON. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If an argument is invalid, as for |
NotImplementedError
|
If this system cannot be saved as JSON (a model that is not a
surpyval parametric one), as a shard carries it: run it with
|
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "t")],
... {"a": {"reliability": E([0.1]), "repairability": E([1.0])}},
... )
>>> len(rbd.shards(100.0, 10_000, seed=1, size=4096))
3
availability_from_chunks ¶
availability_from_chunks(chunks, mc_samples=None, *, allow_gaps=False, control_variate=None, conditional=None)
The result of the simulations of chunks (see
simulate_chunk, and the partials run_shard gives back from
shards), as availability gives it: merged, in order of
their positions in the run, whatever order they come in.
Chunks of simulations 0 to N - 1 give the result of
availability(t_simulation, mc_samples=N, seed=seed, ...): the
same per-simulation values (uptimes, the costs' samples) and
timeline, and the same totals, to the last bit: every total is kept
exactly and rounded once, so it does not depend on how the run was
cut.
The chunks must hold simulations 0 to N - 1 with none
missing, so a lost chunk (a shard that never came back) is not
taken for a smaller run; with allow_gaps=True they give
the result of whichever simulations they hold. With costs, the
result's cost is the simulated cost distribution, as cost
gives it.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
chunks
|
SimulationChunk, or an iterable of them
|
The chunks, their |
required |
mc_samples
|
int
|
The run's number of simulations: the chunks must then hold
simulations |
None
|
allow_gaps
|
bool
|
Take chunks with simulations missing between or before them, for the result of those they hold. By default False. |
False
|
control_variate
|
bool
|
None (the default): the means |
None
|
conditional
|
bool
|
None (the default) or False, as for |
None
|
Returns:
| Type | Description |
|---|---|
AvailabilityResult
|
The result of their simulations. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If the chunks are of different runs, of another system (or of
it saved by another RePyability version), or their simulations
overlap or interleave; or, unless |
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD, SimulationChunk
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "t")],
... {"a": {"reliability": E([0.1]), "repairability": E([1.0])}},
... )
>>> saved = rbd.simulate_chunk(100.0, 0, 200, seed=3).to_json()
>>> chunk = SimulationChunk.from_json(saved)
>>> rbd.availability_from_chunks(chunk).n_simulations
200
compare ¶
compare(other, t_simulation, mc_samples=None, seed=None, *, quantity='availability', confidence=0.95, n_jobs=None, engine='auto', state=None, control_variate=None)
How much better (or worse) this system is than other: exactly
where the exact methods work out both systems' expected values, and
otherwise by simulation with common random numbers.
By default, as availability and cost take their means,
where the exact methods work out both systems' expected
values over the window (mission_availability, and
expected_cost their costs), the difference is theirs, exact,
with no error (method="exact"), and nothing is simulated.
control_variate=False simulates it whatever.
Both systems are simulated mc_samples times over
[0, t_simulation] (every component working and new at the start,
or from state), and in each simulation a
component in the same place in both (the same node name, and the
same names down through nested RBDs) draws the same random numbers
in both: the same failures and repairs where it is modelled the
same way, and matching ones (the same quantiles of its own models)
where it is not. The differences between the two systems' results
then come from how the systems differ rather than from chance, so
their mean is a more precise estimate of the difference than the
difference of two independent simulations of the same size: much
more, when the systems differ in a component's models and share
the rest.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
other
|
RepairableRBD
|
The system to compare with. |
required |
t_simulation
|
float
|
The window each simulation covers. |
required |
mc_samples
|
int
|
The number of simulations of each system, by default 10_000. |
None
|
seed
|
int
|
Seed for a reproducible comparison, by default None: a number
drawn from numpy's global RNG (see |
None
|
quantity
|
str
|
|
'availability'
|
confidence
|
float
|
The confidence level of the interval, by default 0.95. |
0.95
|
n_jobs
|
int
|
Run the simulations in parallel on |
None
|
engine
|
str
|
What runs the simulations: |
'auto'
|
state
|
dict
|
The components' states at the start, as for |
None
|
control_variate
|
bool
|
By default None: the exact difference where the exact methods
work out both systems' expected values. False: the simulated
difference, whatever. (True, which controls |
None
|
Returns:
| Type | Description |
|---|---|
ConfidenceInterval
|
The expected difference (this system's quantity minus
|
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
If a component's draws cannot be replayed from a stream of its
own (a model sampled its own way), or, with
|
ImportError
|
With |
Examples:
Two pumps in parallel, each failing about every 10 hours: how much more of a 100-hour window is the system up if a repair takes 1 hour on average instead of 2?
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def pumps(mttr):
... pump = {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1 / mttr]),
... }
... return RepairableRBD(
... [("s", "p1"), ("s", "p2"), ("p1", "t"), ("p2", "t")],
... {"p1": pump, "p2": pump},
... )
>>> gain = pumps(1.0).compare(pumps(2.0), 100.0)
>>> round(gain.estimate, 4), gain.method
(0.0189, 'exact')
Each pump is up mu / (lambda + mu) + lambda / (lambda + mu) *
exp(-(lambda + mu) t) of the time at t, which the exact
methods average over the window. Simulated with common random
numbers:
>>> simulated = pumps(1.0).compare(
... pumps(2.0), 100.0, mc_samples=2000, seed=0,
... control_variate=False,
... )
>>> round(simulated.estimate, 4), round(simulated.standard_error, 5)
(0.0187, 0.00041)
Two independent runs of 2000 simulations would estimate it with a standard error of about 0.00058.
cost ¶
cost(t_simulation, working_nodes=None, broken_nodes=None, method='p', mc_samples=None, verbose=False, seed=None, *, tolerance=None, confidence=0.95, max_samples=None, antithetic=False, n_jobs=None, engine='auto', state=None, shard_map=None, shard_size=None, control_variate=None, conditional=None)
Simulate the cost of running the system for t_simulation.
Runs the availability simulation with cost accumulation and
returns its CostResult: the distribution
of the window's total cost (samples, mean,
percentile(q)), a confidence interval for the mean
(mean_interval()), and per-category and per-component
breakdowns. In each simulation repair_cost and replace_cost
are charged at every failure of their component (a fresh draw for a
cost distribution), a preventive-maintenance cost at every
preventive replacement, downtime_cost per unit time its
component is down (planned outages included), and
downtime_cost_rate per unit time the system is down.
Only this RBD's own costs count: costs declared inside a nested
RepairableRBD are left out. result.cost_rate converges to
the exact expected_cost_rate as the window grows, which is the
cross-check to reach for.
When nothing is priced but the components' acquisition_cost
(see has_costs), the result has no running cost, exactly, and
the acquisition beside it, without a simulation; when nothing at
all is priced there is no cost model to evaluate, and this
refuses. (The same result is available as availability(...).cost
if you also want the availability outputs from the same
simulations; there it is None when nothing at all is priced.)
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t_simulation
|
float
|
Length of the window each simulation covers, in the time unit of the component models. |
required |
working_nodes
|
Collection[Hashable]
|
Nodes held working for the whole window, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed for the whole window, by default None. |
None
|
method
|
str
|
Evaluate the system state from the minimal path sets ( |
'p'
|
mc_samples
|
int
|
Number of simulations, each giving one sample of the window's total cost, by default 10_000. |
None
|
verbose
|
bool
|
If True, displays a progress bar of the simulations, by default False. |
False
|
seed
|
int
|
Seed for a reproducible run, by default None (see
|
None
|
tolerance
|
float
|
Simulate until the mean cost of a window is known to within
|
None
|
confidence
|
float
|
The confidence level |
0.95
|
max_samples
|
int
|
The most simulations a run to |
None
|
antithetic
|
bool
|
Run the simulations in antithetic pairs (see |
False
|
n_jobs
|
int
|
Run the simulations in parallel on |
None
|
engine
|
str
|
What runs the simulations: |
'auto'
|
state
|
dict
|
Start each simulation from the components' current states
rather than new, as for |
None
|
shard_map
|
callable
|
Run the simulations as shards through this map, wherever it
sends them, as for |
None
|
shard_size
|
int
|
The simulations to a shard, with |
None
|
control_variate
|
bool
|
Control the estimate of the mean cost by the system's exact
twin, as for |
None
|
conditional
|
bool
|
Take the expected cost given the histories of the dependent
modules, as for |
None
|
Returns:
| Type | Description |
|---|---|
CostResult
|
The simulated costs. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If nothing at all is priced (see |
NotImplementedError
|
As for |
ImportError
|
As for |
Examples:
One component with MTTF 10 and MTTR 1, at 100 per repair, and 50 per hour of system downtime:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {
... "c": {
... "reliability": surv.Exponential.from_params([0.1]),
... "repairability": surv.Exponential.from_params([1.0]),
... "repair_cost": 100.0,
... }
... },
... downtime_cost_rate=50.0,
... )
>>> result = rbd.cost(t_simulation=100.0, mc_samples=200, seed=0)
>>> round(result.mean, 2) # expected cost of a 100-hour window: exact
1360.33
>>> round(result.sample_mean, 2) # the 200 windows' own
1387.42
>>> round(result.percentile(90), 2) # 9 windows in 10 cost less
1897.0
>>> round(result.cost_rate, 2), round(rbd.expected_cost_rate(), 2)
(13.6, 13.64)
node_availability ¶
node_availability()
Returns each node's long-run (steady-state) availability.
Exact: MTTF / (MTTF + MTTR) for each component (1.0 if it is
repaired instantly), and a nested RepairableRBD's own
mean_availability; the components alone, as node_sf and
the importance measures list them (the input and output nodes
are no components). These are the availabilities the importance
measures are evaluated at.
A component under age replacement at T renews at a failure or
at a preventive replacement, whichever comes first. By the
renewal-reward theorem its availability is its mean up time per
cycle over the mean cycle length,
integral_0^T R / (integral_0^T R + F(T) * MTTR + R(T) * MTTP),
with MTTP the mean maintenance time (0 for "instant").
Under block replacement every T a unit is replaced at each
multiple of T at which it is up, and in between it fails and is
repaired as usual; a replacement due while it is down is skipped.
The block times at which it is up are its renewals, so its
availability is its mean up time between two of them over their
mean distance apart. That needs the expected up time and number of
failures of an alternating renewal process of lives and repairs
within an interval, and what a repair still going on at a block time
carries into the next: they are computed numerically, with an error
falling as the square of the grid step (about 1e-6 or less). With
instant repair and replacement the unit is always up, and it fails
M(T) times per interval, M the renewal function of its
lives.
A component with hidden failures, a constant failure rate
lambda, inspected every tau with instant tests and instant
repair, is up a fraction (1 - exp(-lambda * tau)) / (lambda *
tau) of the time, about 1 - lambda * tau / 2: it is down, from
a failure until the next inspection, half an interval on average.
Other hidden failures have no exact long-run values here: simulate
them.
With limited repair_crews, a component's is its long-run
probability of being up in the crews' Markov chain (see
mean_availability): waiting for a crew, it is down.
Returns:
| Type | Description |
|---|---|
dict[Hashable, float]
|
Node name -> long-run availability, in |
Raises:
| Type | Description |
|---|---|
NotImplementedError
|
If a component is under block replacement with models its exact
values do not cover (a lifetime that is not a surpyval
parametric model with a density, dead-on-arrival units, repairs
that may never end, or repairs or maintenance far longer than
the interval), or has hidden failures its numerical values do
not cover (see |
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": "instant"},
... },
... )
>>> availability = rbd.node_availability()
>>> {node: round(a, 4) for node, a in availability.items()}
{'a': 0.8333, 'b': 1.0}
system_failure_frequency ¶
system_failure_frequency(working_nodes=None, broken_nodes=None)
Returns the system's long-run failure frequency (failures per unit time), by the Birnbaum/Vesely formula.
In steady state the system failure frequency is
omega_sys = sum_i I_B(i) * omega_i
where I_B(i) is node i's Birnbaum importance evaluated at the
nodes' long-run availabilities (see birnbaum_importance) and
omega_i is node i's long-run failure frequency: 1 / (MTTF_i +
MTTR_i) for a component, F(T) / C for one under age
replacement at T (C its mean renewal cycle, see
expected_cost_rate), its mean failures per renewal cycle over
the cycle's mean length for one under block replacement, and a
nested RepairableRBD's own
system_failure_frequency. Exact for independent repairable
nodes, with no simulation. Each Birnbaum importance is a sum of
products of the node availabilities and unavailabilities, so a small
frequency (a reliable redundant system's, say) keeps its full
relative precision. Every system failure counts, including
the zero-length outages an instantly repaired component causes;
planned outages (preventive maintenance) are not failures.
With limited repair_crews the formula is summed over the states
of the crews' Markov chain (see mean_availability): in each, a
component that is up fails at its constant rate, and takes the
system down if it is critical there.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working (they contribute no failures), by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed (they contribute no failures), by default None. |
None
|
Returns:
| Type | Description |
|---|---|
float
|
Expected number of system failures per unit time, in the long run. |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
If a component is under block replacement with models its exact
values do not cover (see |
Examples:
Two components in series, with MTTFs 5 and 2 and MTTRs of 1:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> round(rbd.system_failure_frequency(), 4) # 2/3 * 1/6 + 5/6 * 1/3
0.3889
>>> round(rbd.system_failure_frequency(working_nodes=["a"]), 4)
0.3333
mean_time_between_failures ¶
mean_time_between_failures(working_nodes=None, broken_nodes=None)
Returns the system's long-run Mean Time Between Failures.
Exact, to full precision however rare the failures:
MTBF = 1 / system_failure_frequency — the mean length of
one full up-down cycle, i.e. MTBF = MUT + MDT (when there are
no planned outages: preventive maintenance that takes time also
ends up periods, so MTBF is then longer). Infinite if the system
never fails (e.g. a node held working keeps it up, or one held
broken keeps it down).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
Returns:
| Type | Description |
|---|---|
float
|
The mean time between system failures, possibly |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
Two components in series, with MTTFs 5 and 2 and MTTRs of 1:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> round(rbd.mean_up_time(), 4) # 1 / (0.2 + 0.5)
1.4286
>>> round(rbd.mean_down_time(), 4)
1.1429
>>> round(rbd.mean_time_between_failures(), 4) # MUT + MDT
2.5714
mean_up_time ¶
mean_up_time(working_nodes=None, broken_nodes=None)
Returns the system's long-run Mean Up Time.
Exact: the mean duration of an uninterrupted working period of the
system (sometimes called the repairable-system MTTF),
MUT = mean_availability / system_failure_frequency, to full
precision however rare the failures. A planned
outage (preventive maintenance that takes time) also ends a working
period, so its frequency is added to the failure frequency. If the
system never goes down it is infinite when the system is up and
0.0 when it is always down.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
Returns:
| Type | Description |
|---|---|
float
|
The mean up time, possibly |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
For exponential failures in series, MUT is one over the sum of the failure rates:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> round(rbd.mean_up_time(), 4) # 1 / (0.2 + 0.5)
1.4286
>>> round(rbd.mean_up_time(working_nodes=["a"]), 4) # b alone
2.0
mean_down_time ¶
mean_down_time(working_nodes=None, broken_nodes=None)
Returns the system's long-run Mean Down Time.
Exact: the mean duration of a system outage (the repairable-system
MTTR), MDT = mean_unavailability / system_failure_frequency,
both worked out in their own right, so a reliable system's keeps
its precision. Planned outages (preventive maintenance that takes
time) count as
outages, so their frequency is added to the failure frequency. If
the system never goes down it is 0.0 when the system is always up
and infinite when it is down.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes always working, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
Returns:
| Type | Description |
|---|---|
float
|
The mean down time, possibly |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
A parallel pair is down only while both are, and the outage ends when either is repaired:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "t"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> round(rbd.mean_down_time(), 4) # 1 / (1.0 + 1.0)
0.5
>>> rbd.mean_down_time(working_nodes=["a"]) # never down
0.0
birnbaum_importance ¶
birnbaum_importance(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None)
Returns the Birnbaum measure of importance for all nodes, evaluated at the nodes' long-run availabilities.
In the guide's Greeks it is delta: how far the system moves with each component (see Sensitivities: the Greeks).
Exact, with no simulation: I_B(i) = A_sys(A_i = 1) -
A_sys(A_i = 0), the system's long-run availability with node i
always working minus that with node i always failed. It is the
probability that node i is critical (the system works if and only
if node i does), and how much the system availability changes per
unit change in node i's availability. The node availabilities are
those of node_availability, with working_nodes and
broken_nodes held at 1 and 0. It is worked out for every node
at once, as the derivative of the system availability: a sum of
products of the node availabilities and unavailabilities (each
worked out in its own right, see mean_unavailability), so a
small one keeps its full relative precision, where the difference
above would cancel in a reliable system.
Note: Birnbaum's measure of importance assumes all nodes are
independent. With common-cause groups, a member's values with it
working and failed are the system's availability given its state
at each long-run time (the groups' Markov chains give the members'
joint states), averaged over the times as every node's are; a node
outside the groups is held up and down, as without them. The other
measures are built on the same values, and beta = 0 gives them
as without the group.
With limited repair_crews, under which a component can wait for
a crew, the components are not independent: node i is then held
working and failed in the crews' Markov chain, which is solved
without it (held either way, it needs no crew), and I_B(i) is
the difference of the system's long-run availabilities so held. The
improvement potential, RAW and RRW are built on the same values;
the criticality and Fussell-Vesely measures are probabilities over
the chain's states. With a crew for each component they are the
independent ones.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes being always available, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
The length of a window |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float or ndarray]
|
Dictionary with node names as keys and Birnbaum importances as
values, for every node except the input and output nodes.
Floats in the long run, at a single time or over a single
window; else arrays in the shape of |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
In series, a node is critical exactly when the others work:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> importance = rbd.birnbaum_importance()
>>> {node: round(i, 4) for node, i in importance.items()}
{'a': 0.6667, 'b': 0.8333}
>>> round(rbd.birnbaum_importance(working_nodes=["a"])["b"], 4)
1.0
improvement_potential ¶
improvement_potential(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None)
Returns the improvement potential of all nodes, evaluated at the nodes' long-run availabilities.
Exact, with no simulation: A_sys(A_i = 1) - A_sys, the gain in
the system's long-run availability if node i never failed, worked
out as I_B(i) * (1 - A_i) (see birnbaum_importance) so that
a small one keeps its precision. The node availabilities are those
of node_availability, with working_nodes and
broken_nodes held at 1 and 0.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes being always available, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
The length of a window |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float or ndarray]
|
Dictionary with node names as keys and improvement potentials as
values, for every node except the input and output nodes.
Floats in the long run, at a single time or over a single
window; else arrays in the shape of |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> potential = rbd.improvement_potential()
>>> {node: round(p, 4) for node, p in potential.items()}
{'a': 0.1111, 'b': 0.2778}
risk_achievement_worth ¶
risk_achievement_worth(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None)
Returns the RAW importance per Modarres & Kaminskiy, evaluated at the nodes' long-run availabilities. That is RAW_i = (unavailability of system given i failed) / (nominal system unavailability).
Exact, with no simulation: RAW_i = U_sys(A_i = 0) / U_sys,
where U = 1 - A is long-run unavailability: the factor by which
the system's unavailability grows if node i is always failed. The
node availabilities are those of node_availability, with
working_nodes and broken_nodes held at 1 and 0. Where the
nominal unavailability is 0 (e.g. every component is repaired
instantly), the ratio is inf or nan, with a numpy warning.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes being always available, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
The length of a window |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float or ndarray]
|
Dictionary with node names as keys and RAW importances as
values, for every node except the input and output nodes.
Floats in the long run, at a single time or over a single
window; else arrays in the shape of |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
In series the system is down whenever a node is, so failing any
node multiplies the unavailability by 1 / U_sys:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> raw = rbd.risk_achievement_worth()
>>> {node: round(r, 4) for node, r in raw.items()}
{'a': 2.25, 'b': 2.25}
risk_reduction_worth ¶
risk_reduction_worth(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None)
Returns the RRW importance per Modarres & Kaminskiy, evaluated at the nodes' long-run availabilities. That is RRW_i = (nominal unavailability of system) / (unavailability of system given i is working).
Exact, with no simulation: RRW_i = U_sys / U_sys(A_i = 1),
where U = 1 - A is long-run unavailability: the factor by which
the system's unavailability would fall if node i never failed. The
node availabilities are those of node_availability, with
working_nodes and broken_nodes held at 1 and 0. It is
inf (with a numpy warning) where making node i perfect would
make the system perfect, as for either node of a parallel pair, and
nan if the system is already never down.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes being always available, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
The length of a window |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float or ndarray]
|
Dictionary with node names as keys and RRW importances as
values, for every node except the input and output nodes.
Floats in the long run, at a single time or over a single
window; else arrays in the shape of |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output node, or is in both sets. |
NotImplementedError
|
As for |
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> rrw = rbd.risk_reduction_worth()
>>> {node: round(r, 4) for node, r in rrw.items()}
{'a': 1.3333, 'b': 2.6667}
criticality_importance ¶
criticality_importance(working_nodes=None, broken_nodes=None, kind='failure', *, x=None, window=None, state=None)
Returns the criticality importance of all nodes, evaluated at the nodes' long-run availabilities.
Exact, with no simulation. With I_B(i) the Birnbaum importance,
A_i node i's availability and A_sys the system's:
kind="failure"(the default) gives the failure-oriented form (Rausand & Høyland),I_B(i) * (1 - A_i) / (1 - A_sys): the probability that node i is down and critical, given that the system is down -- the share of the system's downtime that node i accounts for. It ranks nodes in series by their unavailability. It is computed from the node unavailabilities, through the minimal cut sets, so the system unavailability is not lost to cancellation in1 - A_syshowever available the system is. It isnanif the system is never down.kind="success"gives the success-oriented form,I_B(i) * A_i / A_sys: the probability that node i is working and critical, given that the system is working. It is 1 for every node in series with the rest of the system, so it cannot rank them. It isnanif the system is never up.
The node availabilities are those of node_availability, with
working_nodes and broken_nodes held at 1 and 0.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes being always available, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
kind
|
str
|
|
'failure'
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
The length of a window |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float or ndarray]
|
Dictionary with node names as keys and criticality importances
as values, for every node except the input and output nodes.
Floats in the long run, at a single time or over a single
window; else arrays in the shape of |
Raises:
| Type | Description |
|---|---|
ValueError
|
If a working/broken node is unknown, is the input or output
node, or is in both sets; or if |
NotImplementedError
|
As for |
References
M. Rausand and A. Høyland, System Reliability Theory: Models, Statistical Methods, and Applications, 2nd edition, Wiley, 2004.
Examples:
Two nodes in series, b down twice as often as a: b
accounts for more of the system's downtime (the shares add to less
than 1, as neither is critical while both are down):
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> criticality = rbd.criticality_importance()
>>> {node: round(c, 4) for node, c in criticality.items()}
{'a': 0.25, 'b': 0.625}
The success-oriented form cannot tell them apart:
>>> criticality = rbd.criticality_importance(kind="success")
>>> {node: round(c, 4) for node, c in criticality.items()}
{'a': 1.0, 'b': 1.0}
fussell_vesely ¶
fussell_vesely(fv_type='c', working_nodes=None, broken_nodes=None, method='exact', *, x=None, window=None, state=None)
Calculate Fussell-Vesely importance of all nodes, evaluated at the nodes' long-run availabilities.
The Fussell-Vesely importance of node i is the probability that
some minimal cut set containing node i has failed (all its nodes
down), over the probability that the system has failed: the share
of the system's unavailability that involves node i. Here a node's
probability of having failed is its long-run unavailability,
1 - A_i, from node_availability (with working_nodes and
broken_nodes held at availability 1 and 0), and the system's is
1 - A_sys; both are exact, with no simulation, and worked out
in their own right (see mean_unavailability), so small ones
keep their precision. method="exact" (the default) works out
the probability that some cut set containing node i has failed
exactly, from the exact engine, so the measure is between 0 and 1.
method="rare_event" sums the cut sets' probabilities instead,
the usual rare-event approximation, which with large
unavailabilities can exceed 1. If the system never fails, the ratio
is nan or inf, with a numpy warning.
Typically this measure is implemented using cut-sets as mentioned
above, although it can be implemented using path-sets. Both are
implemented here, selected by fv_type: "c" takes the
minimal cut sets containing node i, "p" the minimal path sets
containing it. Either way a set has failed when all of its nodes
have, and the probability that one of them has (or, with
"rare_event", the sum of their probabilities) is divided by the
system's unavailability.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
fv_type
|
str
|
Dictates the method of calculation, "c" = cut-set and "p" = path-set, by default "c". |
'c'
|
working_nodes
|
Collection[Hashable]
|
Condition on these nodes being always available, by default None. |
None
|
broken_nodes
|
Collection[Hashable]
|
Condition on these nodes being failed, by default None. |
None
|
method
|
str
|
|
'exact'
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
The length of a window |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict[Any, float or ndarray]
|
Dictionary with node names as keys and Fussell-Vesely importances
as values, for every node except the input and output nodes.
Floats in the long run, at a single time or over a single
window; else arrays in the shape of |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
NotImplementedError
|
As for |
Examples:
In series each node is a cut set on its own, so its importance is its unavailability over the system's:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> fv = rbd.fussell_vesely()
>>> {node: round(v, 4) for node, v in fv.items()}
{'a': 0.375, 'b': 0.75}
parameter_sensitivity ¶
parameter_sensitivity(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None, rel_step=None, of='availability', unit_costs=None)
How the system's availability (or its cost rate) moves with each lever: the derivative in each of its components' parameters, and the change one more standby unit or repair crew makes.
In the guide's Greeks it is the levers' deltas (see Sensitivities: the Greeks).
The levers are each component's life and repair models'
parameters ("reliability.alpha", "repairability.beta",
... with surpyval's names), its preventive maintenance's
("preventive.interval", "preventive.threshold",
"preventive.opportunity" and "preventive.duration.<name>"),
its tests' ("inspection.interval", "inspection.coverage",
"inspection.offset" and "inspection.duration.<name>"), its
standby group's ("standby.dormancy_factor",
"standby.switching_probability", and "standby.units", one
more unit), its imperfect repair's ("repair.q"), a common-cause
group's (its members' parameters, moved together, and its model's
"ccf_beta", "ccf_gamma", ...), and the repair crews
("repair_crews", one more, under the key None). levers()
lists them, with their values and ranges, and with_levers
builds the RBD with them moved, as they are moved here.
A continuous lever's derivative is a central difference of the
system's own value with the lever moved by rel_step of its
value either way, the diagram rebuilt with the changed spec:
one-sided where one side is not a valid value (a coverage past 1,
say), NaN where neither is. In the long run (the default) that is
a difference of the exact long-run values, worked out as
unavailabilities so that a small one keeps its precision. At times
x (from new, or from state) or over [0, window), with
independent components, a component's lever moves its own curve
alone, and the system's availability is linear in each
component's: the derivative is the component's Birnbaum importance
at each time (see birnbaum_importance) times the difference of
its own point availability, worked out on the component alone,
and over a window the same averaged on the window's quadrature
points (but for a lever that moves its curve's breaks, an interval
or offset, differenced on the system). With limited repair crews
or common-cause groups the system itself is differenced.
In the long run, a block replacement's or a test's interval of a component whose schedule shares a calendar with others' (other block replacements or tests) would move it off their common calendar, where the long-run values jump (they are averages over the schedules' common period, and the components' outages fall together or apart): its derivative takes the component's schedule as apart from the others', its long-run Birnbaum importance times the change in its own long-run availability (and its own cost rate). A test interval moves the full tests with it, every so many tests as before. A common-cause group's interval there is NaN.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Nodes held working, as for |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed, likewise. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float
|
Take the sensitivity of the mean availability over
|
None
|
state
|
dict or str
|
With |
None
|
rel_step
|
float
|
The step of a continuous lever's difference, relative to its
value (absolute where the value is 0), by default |
None
|
of
|
str or tuple of str
|
|
'availability'
|
unit_costs
|
dict
|
The cost of a unit change of each lever to rank, by |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
|
Raises:
| Type | Description |
|---|---|
ValueError
|
For both |
NotImplementedError
|
As the value differenced refuses ( |
Examples:
A component that fails at the rate 0.01 and is repaired at the
rate 0.1 is up 0.1 / 0.11 of the time in the long run: its
availability falls by 0.1 / 0.11 ** 2 per unit of the failure
rate, about 8.264, and rises by 0.01 / 0.11 ** 2, about
0.8264, per unit of the repair rate:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": {"reliability": E([0.01]), "repairability": E([0.1])}},
... )
>>> sensitivity = rbd.parameter_sensitivity()["c"]
>>> round(sensitivity["reliability.failure_rate"], 3)
-8.264
>>> round(sensitivity["repairability.failure_rate"], 4)
0.8264
levers ¶
levers()
The levers that parameter_sensitivity moves, in the order it
reports them: each component's life and repair models'
parameters, its maintenance's and tests' options, its standby
group's and its imperfect repair's, each common-cause group's, and
the repair crews.
Each is a Lever: whose it is and its name, as
parameter_sensitivity keys them, its value, the range of its
values, whether it is discrete (one more standby unit or repair
crew), and whether it moves a calendar its component shares with
others. Nothing is worked out, so a report can name the levers and
show their values without the sensitivities;
with_levers sets them.
The levers are the components' own models' parameters and
options: those of the units inside a composite life model (a
StandbyModel, a LoadSharingModel, a DegradingNode) or
of a nested diagram are not reached (see the guide's Greeks).
Returns:
| Type | Description |
|---|---|
list of Lever
|
The levers, in |
Examples:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": {"reliability": E([0.01]), "repairability": E([0.1])}},
... )
>>> for lever in rbd.levers():
... print(lever.key, lever.name, lever.value, lever.bounds)
c reliability.failure_rate 0.01 (0.0, inf)
c repairability.failure_rate 0.1 (0.0, inf)
with_levers ¶
with_levers(values)
A copy of this RBD with levers at new values, each moved
as parameter_sensitivity moves it, so that a what-if agrees
with the sensitivities: a test interval takes the full tests with it
(every so many tests, as before), and a common-cause group's
members' parameter moves for every member. levers() lists the
levers.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
values
|
dict
|
Each lever's new value, by its |
required |
Returns:
| Type | Description |
|---|---|
RepairableRBD
|
The new RBD, built as this one was with the levers moved; this one is unchanged. |
Raises:
| Type | Description |
|---|---|
ValueError
|
For a lever the RBD does not have (with the closest name), a
discrete lever given a fraction, or a value the RBD refuses
(one outside the lever's |
TypeError
|
If |
Examples:
Repairs in 8 rather than 10, on average: the availability rises
from 100 / 110 to 100 / 108:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "c"), ("c", "t")],
... {"c": {"reliability": E([0.01]), "repairability": E([0.1])}},
... )
>>> repair = ("c", "repairability.failure_rate")
>>> faster = rbd.with_levers({repair: 1 / 8})
>>> round(float(faster.mean_availability()), 4)
0.9259
>>> float(faster.repairability["c"].mean())
8.0
mean_availability_uncertainty ¶
mean_availability_uncertainty(uncertainty=None, *, n_draws=1000, seed=None, sampling='random')
The long-run availability over plausible models of the components (parameter uncertainty).
A component's models are estimated from data, so their parameters
are uncertain: epistemic uncertainty, about what the models are,
as opposed to the aleatory variability they describe. Each of the
n_draws draws gives every uncertain model a plausible value, and
the diagram, rebuilt with them, its long-run availability, worked
out as mean_availability works it out (exactly, or numerically
with maintenance): the spread of the draws, and its percentiles
(interval), say how well the availability is known.
A component has several models: its life ("reliability"), its
repair ("repairability"), and the durations of its preventive
maintenance and of its tests ("preventive.duration",
"inspection.duration"), named as parameter_sensitivity's
levers. uncertainty maps an input to its uncertainty: a node; a
tuple of nodes of one population, sharing a model, which each draw
gives them alike; or a common-cause group (the
CCFGroup itself), whose model's
parameters, or alternative models, are drawn as for
NonRepairableRBD.sf_uncertainty. A node may come under several
inputs, a role of it under one: a fleet's life fit, shared by two
pumps whose repairs were recorded apart, is {("a", "b"):
{"reliability": "fit"}, "a": {"repairability": "fit"}, "b":
{"repairability": "fit"}}, the life drawn once a draw. An
input drawn without other nodes that hold the same model object,
which then keep it as fitted, is warned about. A node's uncertainty
is
"fit": every one of its models that is a surpyval fit with a parameter covariance (covariance()), drawn from its normal approximation (on the log scale for a positive parameter, the logit scale for one in (0, 1));{role: uncertainty}: each model named drawn as its uncertainty says ("fit",{parameter: distribution}or a list of models), the others kept;{parameter: distribution}or a list of models: its life's.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
uncertainty
|
dict
|
|
None
|
n_draws
|
int
|
The number of draws, by default 1000: one evaluation each. |
1000
|
seed
|
int
|
Seed for the draws. |
None
|
sampling
|
str
|
|
'random'
|
Returns:
| Type | Description |
|---|---|
UncertaintyResult
|
The availability of every draw ( |
Raises:
| Type | Description |
|---|---|
ValueError
|
If |
Examples:
A pump whose repair time was fitted to 12 repairs, in series with a valve: the repair's uncertainty is the availability's.
>>> import numpy as np
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> repairs = surv.LogNormal.fit(
... np.exp(1.5 + 0.4 * np.linspace(-1.6, 1.6, 12))
... )
>>> rbd = RepairableRBD(
... [("s", "pump"), ("pump", "valve"), ("valve", "t")],
... {
... "pump": {
... "reliability": surv.Exponential.from_params([0.01]),
... "repairability": repairs,
... },
... "valve": {
... "reliability": surv.Exponential.from_params([0.002]),
... "repairability": surv.Exponential.from_params([0.5]),
... },
... },
... )
>>> result = rbd.mean_availability_uncertainty(n_draws=2000, seed=0)
>>> round(result.nominal, 4)
0.9499
>>> lower, upper = result.interval(0.9)
>>> round(lower, 4), round(upper, 4)
(0.9392, 0.9577)
point_availability_uncertainty ¶
point_availability_uncertainty(x, uncertainty=None, *, n_draws=1000, seed=None, sampling='random', state=None)
The point availability at times x (from new, or from
state) over plausible models of the components, each draw
worked out as point_availability works it out: the inputs and
their uncertainty as for mean_availability_uncertainty.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array - like
|
Times. |
required |
uncertainty
|
dict
|
As for |
None
|
n_draws
|
int
|
The number of draws, by default 1000: one evaluation of the availability over time each. |
1000
|
seed
|
int
|
Seed for the draws. |
None
|
sampling
|
str
|
|
'random'
|
state
|
dict or str
|
The components' current states, as for |
None
|
Returns:
| Type | Description |
|---|---|
UncertaintyResult
|
Per draw, the availability at each time ( |
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
mission_availability_uncertainty ¶
mission_availability_uncertainty(t, uncertainty=None, *, n_draws=1000, seed=None, sampling='random', state=None)
The mission availability over [0, t] (from new, or from
state) over plausible models of the components, each draw
worked out as mission_availability works it out: the inputs and
their uncertainty as for mean_availability_uncertainty.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
t
|
float or array - like
|
The missions' lengths. |
required |
uncertainty
|
dict
|
As for |
None
|
n_draws
|
int
|
The number of draws, by default 1000. |
1000
|
seed
|
int
|
Seed for the draws. |
None
|
sampling
|
str
|
|
'random'
|
state
|
dict or str
|
The components' current states, as for
|
None
|
Returns:
| Type | Description |
|---|---|
UncertaintyResult
|
Per draw, the mission availability of each mission, the
|
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
expected_cost_rate_uncertainty ¶
expected_cost_rate_uncertainty(uncertainty=None, *, n_draws=1000, seed=None, sampling='random')
The long-run cost rate over plausible models of the components,
each draw worked out as expected_cost_rate works it
out: the inputs and their uncertainty as for
mean_availability_uncertainty.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
uncertainty
|
dict
|
As for |
None
|
n_draws
|
int
|
The number of draws, by default 1000. |
1000
|
seed
|
int
|
Seed for the draws. |
None
|
sampling
|
str
|
|
'random'
|
Returns:
| Type | Description |
|---|---|
UncertaintyResult
|
The cost rate of every draw, the |
Raises:
| Type | Description |
|---|---|
(ValueError, NotImplementedError)
|
As for |
uncertainty_importance ¶
uncertainty_importance(x=None, uncertainty=None, *, of='mean_availability', method='delta', n_draws=1000, seed=None, sampling='random', rel_step=None, state=None)
Which input's parameter uncertainty makes the availability (or
the cost rate) uncertain: each uncertain input's share of its
variance (vega), the inputs as for
mean_availability_uncertainty.
In the guide's Greeks it is vega: whose uncertainty widens the answer (see Sensitivities: the Greeks).
of is the quantity: "mean_availability" (the default, the
long run), "point_availability" at the times x,
"mission_availability" over [0, x], or
"expected_cost_rate" (or "cost_rate", as
parameter_sensitivity names it).
method="delta" (the default) linearises the quantity in the
inputs' parameters: Var(Q) ~ sum_k g_k^T Sigma_k g_k, g_k
its derivatives in input k's parameters (parameter_
sensitivity's, of every one of its uncertain models; a
population's summed over its nodes, which move together) and
Sigma_k their covariance (a fit's covariance(), or the
variances of the distributions given). The inputs are independent,
so each one's part is its own term, and the shares add up to 1. A
list of models has no parameters to move: it needs
method="sobol".
method="sobol" draws the inputs as
mean_availability_uncertainty does and estimates the first-order
and total Sobol indices: n_draws draws of two independent sets,
and one more set for each input with it taken from the second
(Jansen's estimators), n_draws * (inputs + 2) evaluations in
all, at the cost of sampling error, which sampling="sobol"
shrinks.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array - like
|
The times for |
None
|
uncertainty
|
dict
|
As for |
None
|
of
|
str
|
The quantity (see above). |
'mean_availability'
|
method
|
str
|
|
'delta'
|
n_draws
|
int
|
With |
1000
|
seed
|
int
|
With |
None
|
sampling
|
str
|
With |
'random'
|
rel_step
|
float
|
With |
None
|
state
|
dict or str
|
For the values over time, the components' current states. |
None
|
Returns:
| Type | Description |
|---|---|
UncertaintyImportance
|
The |
Raises:
| Type | Description |
|---|---|
ValueError
|
For an unknown |
Examples:
A pump whose repair was fitted to 12 repairs, in series with a valve failing half as often, whose repair was fitted to 8: the pump's, down twice as often, makes three quarters of the availability's uncertainty.
>>> import numpy as np
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def fitted(median, n):
... spread = 0.4 * np.linspace(-1.6, 1.6, n)
... return surv.LogNormal.fit(median * np.exp(spread))
>>> rbd = RepairableRBD(
... [("s", "pump"), ("pump", "valve"), ("valve", "t")],
... {
... "pump": {
... "reliability": surv.Exponential.from_params([0.01]),
... "repairability": fitted(4.5, 12),
... },
... "valve": {
... "reliability": surv.Exponential.from_params([0.005]),
... "repairability": fitted(4.0, 8),
... },
... },
... )
>>> parts = rbd.uncertainty_importance()
>>> {n: round(s, 2) for n, s in parts.first_order.items()}
{'pump': 0.74, 'valve': 0.26}
differential_importance ¶
differential_importance(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None, over='components', change='uniform', kind='failure', improving=False, groups=None, rel_step=None)
Each node's (or lever's) share of the change in the system's
availability when they all change together: the differential
importance measure (DIM, Borgonovo & Apostolakis, 2001), in
the long run, at times x or over a window.
In the guide's Greeks it is DIM, the shares of a change (see Sensitivities: the Greeks).
DIM_i = dA/dtheta_i dtheta_i / sum_j dA/dtheta_j dtheta_j, so
the shares add up to 1, and a group's share is the sum of its
members' (groups): what share of a possible gain lies in the
repair times, say, or in the maintenance intervals.
change="uniform" moves every theta by as much;
"proportional" each by the same fraction of itself.
Over the nodes (over="components"), theta is each node's
unavailability (kind="failure") or availability
(kind="success"). A uniform change shares out the Birnbaum
importance, either way; a proportional one the criticality
importance of that kind (the failure-oriented one is the
improvement potential over the system's unavailability). Each is
as birnbaum_importance and criticality_importance work it
out, with x, window and state. A node held working or
failed takes no part (its share is 0).
Over the levers (over="parameters"), the derivatives are
parameter_sensitivity's, keyed (key, lever), and a
proportional change moves each lever by the same fraction of its
value. One more standby unit or repair crew is no derivative, and
takes no part.
Effects that oppose (a failure rate that lowers the
availability, a repair rate that raises it) give shares of either
sign, and some beyond 1; where they cancel (an exponential unit's
failure and repair rates, moved in proportion), there is nothing
to share out, and the shares are NaN. improving=True moves
each lever instead the way that raises the availability (a
failure rate down, a repair rate up), so that every share is of a
gain, between 0 and 1: what share of the gain from improving
every lever by the same fraction lies in each. A uniform change
adds the same amount to levers of different units (an interval in
hours, a rate per hour): over levers, the proportional change is
usually the one to ask for.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Nodes held working, as for |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed, likewise. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float or array - like
|
Over |
None
|
state
|
dict or str
|
With |
None
|
over
|
str
|
|
'components'
|
change
|
str
|
|
'uniform'
|
kind
|
str
|
Over the nodes, with a proportional change, which of each
node's probabilities moves in proportion: |
'failure'
|
improving
|
bool
|
Move each lever the way that raises the availability, so that the shares are of a gain (each one's size, shared out), instead of every one up together (the default, the measure as defined). The nodes' shares are the same either way. |
False
|
groups
|
dict
|
|
None
|
rel_step
|
float
|
The levers' step, as for |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
The shares, by node, |
Raises:
| Type | Description |
|---|---|
ValueError
|
For an unknown |
NotImplementedError
|
As they do. |
Examples:
Two units in series, up 10/11 and 4/5 of the time: a uniform change
of their unavailabilities moves the system's by 0.8 dU_a +
(10/11) dU_b (the Birnbaum importances):
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.1]), "repairability": E([1.0])},
... "b": {"reliability": E([0.25]), "repairability": E([1.0])},
... },
... )
>>> shares = rbd.differential_importance()
>>> {node: round(share, 4) for node, share in shares.items()}
{'a': 0.4681, 'b': 0.5319}
A proportional change weighs each by its unavailability, 0.8 /
11 against (10/11) / 5:
>>> shares = rbd.differential_importance(change="proportional")
>>> round(shares["b"], 4)
0.7143
availability_rate ¶
availability_rate(x, working_nodes=None, broken_nodes=None, *, state=None)
How fast the system's availability is changing at each time
x from new (or from the components' state), and which
components are moving it.
In the guide's Greeks it is theta: what is moving the system now (see Sensitivities: the Greeks).
The components failing and recovering independently, the system's point availability is multilinear in theirs, so
dA/dt = sum_i I_B^i(t) dA_i/dt
exactly, I_B^i(t) component i's Birnbaum importance at the
components' point availabilities then: each term is what that
component is doing to the system, pulling it down (negative) as it
wears or fails, or up as it is restored, and the terms add up to
the system's rate. A component's dA_i/dt is the rate of change
of its point availability (see point_availability), by
differences on the grid that is worked out on, so the rates are
numerical: to about 1e-5 of their size, less at 0 where a
life's density is not smooth there. The down times kept off the
grid (a repair or maintenance that starts at a known time) are
differentiated on their own scale, however short, so the
rates just after a scheduled event are as close.
At a scheduled event a component's availability can jump: a block replacement or test that takes it off line, maintenance due at a fixed age from new. The system's jumps there are reported apart, each split among the components that jump together along the straight path between their values before and after, which, the system being multilinear, gives parts that add up to the jump. The rate at the time of a jump is the rate just after it.
With limited repair crews or common-cause groups the components no
longer fail and recover independently, and the parts come from
the crews' or the groups' Markov chains. The system's rate
is p(t) Q u, p(t) the chain's states' probabilities,
Q its generator and u the system up in each state; each
transition is one component's failure or repair (a crew finishing
a repair and taking the next job belongs to the repair that freed
it), so Q splits into each component's transitions, and each
part, p(t) Q_i u, is exact. A nested RBD, with crews of its
own, takes its part as an independent component does. In a
common-cause group, a member's own cause and its repairs are its
part, and the causes that strike more than one member are the
group's, under the tuple of its members. A hidden group's members
are found at their tests, where the system jumps: as tests change
the members' joint states, which are not independent, a jump with
a test in it is split by the Shapley value of each test and node
that changes then, worked out exactly (for independent nodes alone,
the same as the path).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
x
|
float or array - like
|
Times from new (or from |
required |
working_nodes
|
Collection[Hashable]
|
Nodes held working: their part is 0. |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed, likewise. |
None
|
state
|
dict or str
|
The components' current states, as for |
None
|
Returns:
| Type | Description |
|---|---|
RateBreakdown
|
The system's |
Raises:
| Type | Description |
|---|---|
ValueError
|
For a negative or non-finite time, or as
|
NotImplementedError
|
As |
Examples:
Two pumps in parallel, wearing out, and a valve after them that fails at a constant rate. Early on it is the valve that pulls the system down; by 8, it has settled, and the pumps' wear does:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> def unit(life):
... return {
... "reliability": life,
... "repairability": surv.Exponential.from_params([1.0]),
... }
>>> rbd = RepairableRBD(
... [("s", "p1"), ("s", "p2"), ("p1", "v"), ("p2", "v"),
... ("v", "t")],
... {
... "p1": unit(surv.Weibull.from_params([10, 3])),
... "p2": unit(surv.Weibull.from_params([10, 3])),
... "v": unit(surv.Exponential.from_params([0.02])),
... },
... )
>>> rate = rbd.availability_rate([2.0, 8.0])
>>> rate.node_rate["v"].round(5).tolist()
[-0.0026, -1e-05]
>>> rate.node_rate["p1"].round(5).tolist()
[-3e-05, -0.00156]
barlow_proschan_importance ¶
barlow_proschan_importance(working_nodes=None, broken_nodes=None, *, window=None, state=None)
Each component's share of the system's failures: the probability
that a system failure is caused by the component's (Barlow &
Proschan, 1975), in the long run or over a window from new
(or from the components' state).
In the guide's Greeks it is theta, integrated: who caused the failures (see Sensitivities: the Greeks).
A component's failure fails the system when the component is
critical then, which it is with probability I_B^i(t), so the
system's failures are sum_i integral I_B^i(t) dM_i(t),
M_i(t) the component's expected failures by t: each term is
the system failures that component causes, and its share of the sum
is its Barlow-Proschan importance. In the long run it is its share
of system_failure_frequency; over a window, of
expected_failures (the exact counterpart of the simulated
failure_criticality_index.per_system_failure). Failures of
several components at one instant (exact lifetimes, units dead on
arrival) are split among them as availability_rate splits a
jump.
With common-cause groups, a cause that strikes several members at once is counted for the group, under the tuple of its members, and one that strikes a member alone for that member; with limited repair crews, the shares come from the crews' chain, its components' failures at their rates in each state, over a window from its transient probabilities.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Nodes held working (they cause nothing). |
None
|
broken_nodes
|
Collection[Hashable]
|
Nodes held failed (likewise). |
None
|
window
|
float or array - like
|
Over |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
Each component's share (and each common-cause group's, under the tuple of its members), adding up to 1: floats in the long run or for one window, else arrays in the windows' shape. NaN where the system cannot fail. |
Raises:
| Type | Description |
|---|---|
ValueError
|
For |
NotImplementedError
|
As |
Examples:
Two units in series, failing at rates 0.2 and 0.5 and repaired at rate 1: in the long run a share of the system's failures in proportion to each one's rate, while the system is up:
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> rbd = RepairableRBD(
... [("s", "a"), ("a", "b"), ("b", "t")],
... {
... "a": {"reliability": E([0.2]), "repairability": E([1.0])},
... "b": {"reliability": E([0.5]), "repairability": E([1.0])},
... },
... )
>>> shares = rbd.barlow_proschan_importance()
>>> {node: round(share, 4) for node, share in shares.items()}
{'a': 0.2857, 'b': 0.7143}
joint_importance ¶
joint_importance(working_nodes=None, broken_nodes=None, *, x=None, window=None, state=None)
The joint (second-order) importance of each pair of components:
whether improving the two together is worth more than improving
each, in the long run, at times x or over a window.
In the guide's Greeks it is gamma: complements or substitutes (see Sensitivities: the Greeks).
JRI(i, j) = d2A / dA_i dA_j = A(1_i, 1_j) - A(1_i, 0_j) -
A(0_i, 1_j) + A(0_i, 0_j), the system's availability with
component i up or down and j up or down: how much j's
Birnbaum importance rises when i goes from down to up. Positive,
the two are complements, as in series: improving either makes
improving the other worth more. Negative, they are substitutes, as
in parallel. It is each component's Birnbaum importance with the
other held up, less that with it held down, as
birnbaum_importance works them out, with x, window and
state; with limited repair crews, each pair held in the crews'
chain. With common-cause groups, the system's availability given
each of the pair up and down is conditioned on their states, as
birnbaum_importance conditions on one: two members of one
group on their joint state, from the group's chain, and a
component outside the groups held (see Common-cause
failures).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
working_nodes
|
Collection[Hashable]
|
Components held up: their pairs are 0. |
None
|
broken_nodes
|
Collection[Hashable]
|
Components held down, likewise. |
None
|
x
|
float or array - like
|
Times from new (or from |
None
|
window
|
float
|
Over |
None
|
state
|
dict or str
|
With |
None
|
Returns:
| Type | Description |
|---|---|
dict
|
|
Raises:
| Type | Description |
|---|---|
ValueError
|
As |
NotImplementedError
|
As |
Examples:
Two units in parallel, then a third, each up 10/11 of the time: the pair are substitutes, and each is a complement of the third.
>>> import surpyval as surv
>>> from repyability import RepairableRBD
>>> E = surv.Exponential.from_params
>>> unit = {"reliability": E([0.1]), "repairability": E([1.0])}
>>> rbd = RepairableRBD(
... [("s", "a"), ("s", "b"), ("a", "c"), ("b", "c"), ("c", "t")],
... {"a": unit, "b": unit, "c": unit},
... )
>>> {pair: round(v, 4) for pair, v in rbd.joint_importance().items()}
{('a', 'b'): -0.9091, ('a', 'c'): 0.0909, ('b', 'c'): 0.0909}