Defining Fault Sampling Approaches in fmdtools

Fault Sampling is used to evaluate the overall resilience of a system to a set of faults and the corresponding risks associated with these faults. There is no single best way to define the set of scenarios to evaluate resilience with, because a given resilience analysis may need more or less detail, support more or less computational time, or be interested in specific scenario types of interest.

Thus, there are a number of use-cases supported by fmdtools for different sampling models. This document will demonstrate and showcase a few of them.

Copyright © 2024, United States Government, as represented by the Administrator of the National Aeronautics and Space Administration. All rights reserved.

The “"Fault Model Design tools - fmdtools version 2"” software is licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

from fmdtools.sim.sample import SampleApproach, FaultSample, FaultDomain
from fmdtools.analyze.tabulate import FMEA, Comparison
import fmdtools.sim.propagate as prop
import numpy as np

Basics

Fault sampling involves: - Defining faults and fault models for each function/component/flows of the model is done using the ‘Mode’ Class - Defining a fault sampling approach, using the: - SampleApproach class, or - FaultDomain and FaultSample calsses directly - Propagating faults through the model (using the propagate.faultsample method) Before proceeding, it can be helpful to look through their respective documentation.

Model Setup

Consider the rover model in rover_model

import inspect
from rover_model import Power, PowerMode

This rover has a Power function:

print(inspect.getsource(Power))

class Power(Function):
    """Rover power supply."""

    __slots__ = ("ee_15", "ee_5", "ee_12", "switch")
    container_s = PowerState
    container_m = PowerMode
    flow_ee_15 = EE
    flow_ee_5 = EE
    flow_ee_12 = EE
    flow_switch = Switch

    def static_behavior(self, time):
        """Determine power use based on mode."""
        if self.m.in_mode("off"):
            self.off_power()
        elif self.m.in_mode("supply"):
            self.supply_power()
        elif self.m.in_mode("short"):
            self.short_power()
        elif self.m.in_mode("no_charge"):
            self.no_charge_power()

        if self.m.in_mode("charge"):
            self.charge_power_usage()
        else:
            self.power_usage()
            if self.m.in_mode("short"):
                self.short_power_usage()

    def dynamic_behavior(self, time):
        """Charge increment over time."""
        self.s.inc(charge=-self.s.power / 100)
        self.s.limit(charge=(0, 100))

    def short_power(self):
        """Power in case of a short has normal voltage."""
        self.ee_5.s.v = 5
        self.ee_12.s.v = 12
        self.ee_15.s.v = 15

    def no_charge_power(self):
        """Battery is out of charge."""
        self.ee_5.s.v = 0
        self.ee_12.s.v = 0
        self.ee_15.s.v = 0

    def off_power(self):
        """Power supply is shut off."""
        self.ee_5.s.put(v=0, a=0)
        self.ee_12.s.put(v=0, a=0)
        self.ee_15.s.put(v=0, a=0)
        if self.switch.s.power:
            self.m.set_mode("supply")

    def supply_power(self):
        """Power supply is in supply mode."""
        if self.s.charge > 0:
            self.ee_5.s.v = 5
            self.ee_12.s.v = 12
            self.ee_15.s.v = 15
        else:
            self.m.set_mode("no_charge")
        if not self.switch.s.power:
            self.m.set_mode("off")

    def power_usage(self):
        """Calculate the power usage in general."""
        self.s.power = (1.0 + self.ee_12.s.mul("v", "a") +
                        self.ee_5.s.mul("v", "a") + self.ee_15.s.mul("v", "a"))

    def charge_power_usage(self):
        """Calculate the power usage when the battery charges."""
        self.s.power = -1
        if self.s.charge == 100:
            self.m.set_mode("off")

    def short_power_usage(self):
        """Calculate power usage when there is a short (calculated as double)."""
        self.s.power = self.s.power * 2
        if self.s.charge == 0:
            self.m.set_mode("no_charge")
        if not self.switch.s.power:
            self.m.set_mode("off")

Which contains the mode PowerMode:

print(inspect.getsource(PowerMode))

class PowerMode(Mode):
    """
    Possible modes for Power function.

    Modes
    -------
    no_charge : Fault
        Battery is out of charge.
    short: Fault
        There is a short.
    supply: Mode
        supply power
    charge: Mode
        charge battery
    standby: Mode
        power supply is in stand by
    off: Mode
        power supply is off
    """

    fm_args = {"no_charge": (1e-5, 100, {"off": 1.0}),
               "short": (1e-5, 100, {"supply": 1.0})}
    opermodes = ("supply", "charge", "off")
    mode: str = "off"
    exclusive = True

The class variable fm_args specifies that there are two possible modes to inject, “no_charge”, and “shortn”, with no more information given for each mode. More information has been added in the tuples of the dictionary, including: - rate - repair cost - phase dictionary The phase dictionary is important because it specifies that this mode is to occur in a given phase. In this case, no_charge is supposed to only occur during the standby phase while short is only supposed to occur during the supply phase. In this Mode, these phases correspond to the operational modes (opermodes), but they may correspond to other operational modes also.

All of these fields are optional, but they do help us develop a more informed statistical sample of the fault modes.

Setting up a FaultSample

Sampling using FaultSample first requires creating a FaultDomain to sample from. These faultdomains can be created from both models and individual functions:

fd_power = FaultDomain(Power())
fd_power.add_all()
fd_power

FaultDomain with faults:
 -('power', 'no_charge')
 -('power', 'short')

from rover_model import Rover
fd_rvr = FaultDomain(Rover())
fd_rvr.add_all()
fd_rvr

FaultDomain with faults:
 -('power', 'no_charge')
 -('power', 'short')
 -('perception', 'bad_feed')
 -('plan_path', 'no_con')
 -('plan_path', 'crash')
 -('drive', 'hmode_0')
 -('drive', 'hmode_1')
 -('drive', 'hmode_2')
 -('drive', 'hmode_3')
 -('drive', 'hmode_4')
 -...more

Note that there are several methods in FaultDomain which let us specify the list of faults we want to sample from, e.g.:

fd_short = FaultDomain(Rover())
fd_short.add_all_modes("short")
fd_short

FaultDomain with faults:
 -('power', 'short')

or:

fd_pwr = FaultDomain(Rover())
fd_pwr.add_all_fxn_modes("power")
fd_pwr

FaultDomain with faults:
 -('power', 'no_charge')
 -('power', 'short')

We can then sample this domain using a FaultSample:

fs_pwr = FaultSample(fd_pwr, def_mdl_phasemap=False)

Note that FaultSamples have two main variables: faultdomain and phasemap. A PhaseMap is essentially a dictionary of phases to sample from.

In the above case, we don’t want to use phase information to form the sample, so we don’t provide any and we set def_mdl_phasemap=False, since this would get phase information from the model.

To add scenarios to the FaultSample, we can then use one of the add methods:

fs_pwr.add_fault_times([2,5,10])
fs_pwr

FaultSample of scenarios:
 - power_no_charge_t2
 - power_no_charge_t5
 - power_no_charge_t10
 - power_short_t2
 - power_short_t5
 - power_short_t10

As shown, this adds the list of faults in the faultdomain over the given times.

Note the underlying rate information in these scenarios is all the same:

fs_pwr.scenarios()

[SingleFaultScenario(sequence={2.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(2,), function='power', fault='no_charge', rate=1e-05, name='power_no_charge_t2', time=2, phase=''),
 SingleFaultScenario(sequence={5.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(5,), function='power', fault='no_charge', rate=1e-05, name='power_no_charge_t5', time=5, phase=''),
 SingleFaultScenario(sequence={10.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(10,), function='power', fault='no_charge', rate=1e-05, name='power_no_charge_t10', time=10, phase=''),
 SingleFaultScenario(sequence={2.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(2,), function='power', fault='short', rate=1e-05, name='power_short_t2', time=2, phase=''),
 SingleFaultScenario(sequence={5.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(5,), function='power', fault='short', rate=1e-05, name='power_short_t5', time=5, phase=''),
 SingleFaultScenario(sequence={10.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(10,), function='power', fault='short', rate=1e-05, name='power_short_t10', time=10, phase='')]

But we know better than this–some of the faults should have zero rate if they are going to occur in phases they don’t apply to!

If we want to sample the given phases for the sample, we can additionally pass a phasemap generated by first running the model in the nominal state:

res, hist = prop.nominal(Rover())

We can then get phase information from this history using the fmdtools.analyze.phases.from_hist:

from fmdtools.analyze.phases import from_hist, phaseplot
phases = from_hist(hist)
phases

{'power': PhaseMap({'off': [0.0, 0.0], 'supply': [1.0, 112.0]}, {'off': {'off'}, 'supply': {'supply'}}),
 'perception': PhaseMap({'off': [0.0, 1.0], 'feed': [2.0, 112.0]}, {'off': {'off'}, 'feed': {'feed'}}),
 'plan_path': PhaseMap({'standby': [0.0, 4.0], 'drive': [5.0, 112.0]}, {'standby': {'standby'}, 'drive': {'drive'}}),
 'override': PhaseMap({'off': [0.0, 1.0], 'standby': [2.0, 112.0]}, {'off': {'off'}, 'standby': {'standby'}})}

Which can be visualized using:

fig = phaseplot(phases)

The PhaseMap for Power is here in power:

phases['power']

PhaseMap({'off': [0.0, 0.0], 'supply': [1.0, 112.0]}, {'off': {'off'}, 'supply': {'supply'}})

Which we can use to create a FaultSample which only samples the phases corresponding to the information given in Mode:

fs_pwr = FaultSample(fd_pwr, phasemap=phases['power'])
fs_pwr.add_fault_phases()
fs_pwr

FaultSample of scenarios:
 - power_no_charge_t0p0
 - power_short_t0p0
 - power_no_charge_t56p0
 - power_short_t56p0

If we look at the rate information, however:

fs_pwr.scenarios()

[SingleFaultScenario(sequence={0.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(0.0,), function='power', fault='no_charge', rate=1e-05, name='power_no_charge_t0p0', time=0.0, phase='off'),
 SingleFaultScenario(sequence={0.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(0.0,), function='power', fault='short', rate=0.0, name='power_short_t0p0', time=0.0, phase='off'),
 SingleFaultScenario(sequence={56.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(56.0,), function='power', fault='no_charge', rate=0.0, name='power_no_charge_t56p0', time=56.0, phase='supply'),
 SingleFaultScenario(sequence={56.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(56.0,), function='power', fault='short', rate=1e-05, name='power_short_t56p0', time=56.0, phase='supply')]

The rate for scenarios outside the phases is zero!

We can remove these scenarios using FaultSample.prune_scenarios:

fs_pwr.prune_scenarios("rate", np.greater, 0.0)
fs_pwr

FaultSample of scenarios:
 - power_no_charge_t0p0
 - power_short_t56p0

fs_pwr.scenarios()

[SingleFaultScenario(sequence={0.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(0.0,), function='power', fault='no_charge', rate=1e-05, name='power_no_charge_t0p0', time=0.0, phase='off'),
 SingleFaultScenario(sequence={56.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(56.0,), function='power', fault='short', rate=1e-05, name='power_short_t56p0', time=56.0, phase='supply')]

As shown, now the only scenarios in the FaultSample are ones which have nonzero rate.

assert all([scen.rate > 0 for scen in fs_pwr.scenarios()])

To enable multiple samples to be generated for different faultdomains accross the model, we can use SampleApproach, e.g.:

phases

{'power': PhaseMap({'off': [0.0, 0.0], 'supply': [1.0, 112.0]}, {'off': {'off'}, 'supply': {'supply'}}),
 'perception': PhaseMap({'off': [0.0, 1.0], 'feed': [2.0, 112.0]}, {'off': {'off'}, 'feed': {'feed'}}),
 'plan_path': PhaseMap({'standby': [0.0, 4.0], 'drive': [5.0, 112.0]}, {'standby': {'standby'}, 'drive': {'drive'}}),
 'override': PhaseMap({'off': [0.0, 1.0], 'standby': [2.0, 112.0]}, {'off': {'off'}, 'standby': {'standby'}})}

sa = SampleApproach(Rover(), phasemaps=phases)
# adding fault domains
sa.add_faultdomain("drive", "fault", "drive", "hmode_1")
sa.add_faultdomain("plan_path", "all_fxn_modes", "plan_path")
sa.add_faultdomain("power", "all_fxn_modes", "power")
sa.add_faultsample("drive", "fault_phases", "drive", phasemap="plan_path")
sa.add_faultsample("plan_path", "fault_phases", "plan_path", phasemap="plan_path")
sa.add_faultsample("power", "fault_phases", "power", phasemap="power")
sa.prune_scenarios()

sa

SampleApproach for rover with:
 faultdomains: drive, plan_path, power
 faultsamples: drive, plan_path, power

sa.scenarios()

[SingleFaultScenario(sequence={58.0: Injection(faults={'drive': ['hmode_1']}, disturbances={})}, times=(58.0,), function='drive', fault='hmode_1', rate=0.02857142857142857, name='drive_hmode_1_t58p0', time=58.0, phase='drive'),
 SingleFaultScenario(sequence={2.0: Injection(faults={'plan_path': ['no_con']}, disturbances={})}, times=(2.0,), function='plan_path', fault='no_con', rate=0.0001, name='plan_path_no_con_t2p0', time=2.0, phase='standby'),
 SingleFaultScenario(sequence={2.0: Injection(faults={'plan_path': ['crash']}, disturbances={})}, times=(2.0,), function='plan_path', fault='crash', rate=0.0001, name='plan_path_crash_t2p0', time=2.0, phase='standby'),
 SingleFaultScenario(sequence={58.0: Injection(faults={'plan_path': ['no_con']}, disturbances={})}, times=(58.0,), function='plan_path', fault='no_con', rate=0.0001, name='plan_path_no_con_t58p0', time=58.0, phase='drive'),
 SingleFaultScenario(sequence={58.0: Injection(faults={'plan_path': ['crash']}, disturbances={})}, times=(58.0,), function='plan_path', fault='crash', rate=0.0001, name='plan_path_crash_t58p0', time=58.0, phase='drive'),
 SingleFaultScenario(sequence={0.0: Injection(faults={'power': ['no_charge']}, disturbances={})}, times=(0.0,), function='power', fault='no_charge', rate=1e-05, name='power_no_charge_t0p0', time=0.0, phase='off'),
 SingleFaultScenario(sequence={56.0: Injection(faults={'power': ['short']}, disturbances={})}, times=(56.0,), function='power', fault='short', rate=1e-05, name='power_short_t56p0', time=56.0, phase='supply')]

This is mostly useful when we would like to sample different functions in a model differently than others (e.g., using different phases), but still want to propagate the scenarios together as a part of a single sample.

Propagating Faults

Given the FaultSample approach, faults can then be propagated through the model to get results. Note that these faults can be sampled in parallel if desired using a user-provided pool (see the parallel pool tutorial in the \pump example folder).

res, hist = prop.fault_sample(Rover(), fs_pwr)

SCENARIOS COMPLETE: 100%|██████████| 2/2 [00:00<00:00,  8.89it/s]

res

power_no_charge_t0p0.endclass.rate: 1e-05
power_no_charge_t0p0.endclass.cost:    0
power_no_charge_t0p0.endclass.prob:  1.0
power_no_charge_t0p0                   0
power_no_charge_t0p0.endclass.in_bound: True
power_no_charge_t0p0               False
power_no_charge_t0p0                   1
power_no_charge_t0p0                   1
power_no_charge_t0p0.endclass.end_dist: 28.91724385770989
power_no_charge_t0p0                 0.0
power_no_charge_t0p0.endclass.faults: {'power': ['no_charge']}
power_no_charge_t0p0incomplete mission faulty
power_no_charge_t0p0.endclass.end_x: 0.0
power_no_charge_t0p0.endclass.end_y: 0.0
power_no_charge_t0p0.endclass.endpt: array(2)
power_short_t56p0.endclass.rate:   1e-05
power_short_t56p0.endclass.cost:       0
power_short_t56p0.endclass.prob:     1.0
power_short_t56p0.en                   0
power_short_t56p0.endclass.in_bound: True
power_short_t56p0.endclass.at_finish: True
power_short_t56p0.endclass.line_dist:  1
power_short_t56p0.endclass.num_modes:  1
power_short_t56p0.endclass.end_dist: 0.0
power_short_t56p0.en 0.11859644458085598
power_short_t56p0.endclass.faults: {'power': ['short']}
power_short_t56p0.en              faulty
power_short_t56p0.endclass.end_x: 29.254775331608617
power_short_t56p0.endclass.end_y: -0.7827640334587979
power_short_t56p0.endclass.endpt: array(2)
nominal.endclass.rate:               1.0
nominal.endclass.cost:                 0
nominal.endclass.prob:               1.0
nominal.endclass.expected_cost:        0
nominal.endclass.in_bound:          True
nominal.endclass.at_finish:         True
nominal.endclass.line_dist:            1
nominal.endclass.num_modes:            0
nominal.endclass.end_dist:           0.0
nominal.endclass.tot_deviation: 0.11859644458085598
nominal.endclass.faults:              {}
nominal.endclass.classification: nominal mission
nominal.endclass.end_x: 29.254775331608617
nominal.endclass.end_y: -0.7827640334587979
nominal.endclass.endpt:         array(2)

These responses can be visualized over the given faults:

fmea = FMEA(res, fs_pwr, metrics = ["end_dist", "line_dist", "tot_deviation"])
fmea.as_table()

		end_dist	line_dist	tot_deviation
power	short	0.000000	1	0.118596
	no_charge	28.917244	1	0.000000

fmea.as_plots(cols=2)

(<Figure size 600x400 with 4 Axes>,
 array([<Axes: title={'center': 'end_dist'}>,
        <Axes: title={'center': 'line_dist'}, xlabel="['function', 'fault']">,
        <Axes: title={'center': 'tot_deviation'}, xlabel="['function', 'fault']">,
        <Axes: >], dtype=object))

Or over time/any other variable:

comp = Comparison(res, fs_pwr, metrics = ["end_dist", "line_dist", "tot_deviation"], factors =['time'])
comp.as_table()

	end_dist	line_dist	tot_deviation
56.0	0.000000	0.00001	0.000001
0.0	0.000289	0.00001	0.000000

comp.as_plots(cols=2)

(<Figure size 600x400 with 4 Axes>,
 array([<Axes: title={'center': 'end_dist'}>,
        <Axes: title={'center': 'line_dist'}, xlabel='time'>,
        <Axes: title={'center': 'tot_deviation'}, xlabel='time'>, <Axes: >],
       dtype=object))