Urban Drone Demo
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.
[1]:
import fmdtools.sim.propagate as propagate
from fmdtools.analyze import tabulate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
import numpy as np
Model Overview
The drone model is defined in drone_mdl_urban.py, along with some visualization functions.
[2]:
from drone_mdl_urban import Drone
[3]:
mdl = Drone()
This is the model structure:
[4]:
mg = mdl.as_modelgraph()
fig, ax = mg.draw()
We can also view the grid environment using its show methods:
[5]:
collections={"all_occupied": {"color": "red"}, "all_allowed": {"color": "blue"}, "start": {"color": "blue"}, "end": {"color": "blue"}}
fig, ax = mdl.flows['environment'].c.show({"height": {}}, collections=collections)
Which shows the Start, End, and allowed/unsafe locations in the 1000x1000-m grid. In this display, line thickness corresponds to building height, and hatching corresponds to whether or not the space is occupied. We can also display this using show.coord3d:
[6]:
fig, ax = mdl.flows['environment'].c.show_z("height", collections=collections)
Nominal Simulation
Below we show how this drone performs in the nominal scenario.
[7]:
results_nom, hist_nom =propagate.nominal(mdl)
fig, axs = hist_nom.plot_line("fxns.store_ee.s.soc",
'flows.ee_1.s.rate',
'flows.dofs.s.z',
'fxns.plan_path.m.mode',
'fxns.plan_path.s.pt',
'flows.des_traj.s.power')
As shown, the flight ends fairly quickly (in 10 minutes), with the drone successively proceeding through points in the flight plan.
We can also view this flightpath in 3-d space using History.plot_trajectories:
[8]:
fig, ax = hist_nom.plot_trajectories('dofs.s.x', 'dofs.s.y', 'dofs.s.z')
Trajectory plots can be overlaid on top of environment plots. In this case, we defined the method plot_env_with_traj and plot_env_with_traj_z for this case.
[9]:
from drone_mdl_urban import plot_env_with_traj, plot_env_with_traj_z
[10]:
plot_env_with_traj_z(hist_nom, mdl)
[10]:
(<Figure size 400x400 with 1 Axes>,
<Axes3D: title={'center': 'trajectory'}, xlabel='x', ylabel='y'>)
[11]:
plot_env_with_traj(hist_nom, mdl)
[11]:
(<Figure size 500x500 with 2 Axes>,
<Axes: title={'center': 'trajectory'}, xlabel='x', ylabel='y'>)
As shown, this is a rather simple straight-line path. If we wanted a more complex scenario, we could make the path more complex by adding multiple destinations or planning the path based on allowed flight/landing locations.
As it is, we may also want to adjust the timestep/speed to get more resolution, since the drone only has a few discrete timesteps in the air.
The results for the simulation are:
[12]:
results_nom
[12]:
endclass:
--rate: 1.0
--cost: 0.0
--expected_cost: 0.0
--repcost: 0
--unsafe_flight_time: 0
--body_strikes: 0.0
--head_strikes: 0.0
--property_restrictions: 0
--safecost: 0.0
--landcost: 0
--p_safety: 0.0
--severities: {'hazardous': 0.0, 'minor': 1.0}
Resilience model
A number of different faults have been implemented in the system.
For example, here we inject a mechanical fault in the left-rear rotor during flight:
[13]:
mdl.fxns['affect_dof'].m.faultmodes
[13]:
{'lf_short': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'lf_openc': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'lf_ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'lf_ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'lf_ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases={}, units='sim'),
'lf_mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases={}, units='sim'),
'lf_mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases={}, units='sim'),
'lf_propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases={}, units='sim'),
'lf_propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases={}, units='sim'),
'lf_propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases={}, units='sim'),
'lr_short': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'lr_openc': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'lr_ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'lr_ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'lr_ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases={}, units='sim'),
'lr_mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases={}, units='sim'),
'lr_mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases={}, units='sim'),
'lr_propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases={}, units='sim'),
'lr_propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases={}, units='sim'),
'lr_propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases={}, units='sim'),
'rf_short': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'rf_openc': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'rf_ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'rf_ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'rf_ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases={}, units='sim'),
'rf_mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases={}, units='sim'),
'rf_mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases={}, units='sim'),
'rf_propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases={}, units='sim'),
'rf_propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases={}, units='sim'),
'rf_propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases={}, units='sim'),
'rr_short': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'rr_openc': Fault(prob=1.0000000000000002e-06, cost=200, phases={}, units='sim'),
'rr_ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'rr_ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases={}, units='sim'),
'rr_ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases={}, units='sim'),
'rr_mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases={}, units='sim'),
'rr_mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases={}, units='sim'),
'rr_propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases={}, units='sim'),
'rr_propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases={}, units='sim'),
'rr_propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases={}, units='sim')}
[14]:
results_fault, hist_fault =propagate.one_fault(mdl, "affect_dof", "lr_mechbreak", time=2.0)
[15]:
fig, axs = hist_fault.plot_line("fxns.store_ee.s.soc",
'flows.ee_1.s.rate',
'flows.dofs.s.z',
'fxns.plan_path.m.mode',
'fxns.plan_path.s.pt',
'flows.des_traj.s.power')
[16]:
plot_env_with_traj_z(hist_fault, mdl)
[16]:
(<Figure size 400x400 with 1 Axes>,
<Axes3D: title={'center': 'trajectory'}, xlabel='x', ylabel='y'>)
[17]:
plot_env_with_traj(hist_fault, mdl)
[17]:
(<Figure size 500x500 with 2 Axes>,
<Axes: title={'center': 'trajectory'}, xlabel='x', ylabel='y'>)
[18]:
results_fault
[18]:
endclass.rate: 1.0000000000000002e-06
endclass.cost: 28600.0
endclass.expected_cost: 2860.0000000000005
endclass.repcost: 500
endclass.unsafe_flight_time: 281
endclass.body_strikes: 0.0
endclass.head_strikes: 0.0
endclass.property_restrictions: 0
endclass.safecost: 28100.0
endclass.landcost: 0
endclass.p_safety: 0.0
endclass.severities: {'hazardous': 0.0, 'minor': 1.0000000000000002e-06}
Here we inject a large list of faults in the system and evaluate their relative consequences in terms of metrics calculated in find_classification:
[19]:
from fmdtools.analyze.phases import PhaseMap, from_hist
phasemaps = from_hist(hist_nom, fxn_modephases=[])
phasemaps
[19]:
{'ctl_dof': PhaseMap({'nominal': [0.0, 30.0]}, {}),
'plan_path': PhaseMap({'taxi': [0.0, 1.0], 'move': [1.1, 1.6], 'land': [1.7, 2.4], 'taxi1': [2.5, 30.0]}, {})}
[20]:
from drone_mdl_urban import make_move_quad
move_quad=make_move_quad(hist_nom, phasemaps['plan_path'].phases['move'])
move_quad
[20]:
{'samp': 'quadrature',
'quad': {'nodes': [-0.20000000000000018, 1.0], 'weights': [1.0, 0.92]}}
[21]:
from fmdtools.sim.sample import FaultDomain, FaultSample
fd = FaultDomain(mdl)
fd.add_all()
fs = FaultSample(fd, phasemap = phasemaps['plan_path'])
fs.add_fault_phases("move", method = "quad",
args=(move_quad['quad']['nodes'], move_quad['quad']['weights']))
fs
[21]:
FaultSample of scenarios:
- manage_health_lostfunction_t2
- manage_health_lostfunction_t2
- store_ee_nocharge_t2
- store_ee_nocharge_t2
- store_ee_lowcharge_t2
- store_ee_lowcharge_t2
- store_ee_s1p1_short_t2
- store_ee_s1p1_short_t2
- store_ee_s1p1_degr_t2
- store_ee_s1p1_degr_t2
- ... (118 total)
[22]:
endresults, hists = propagate.fault_sample(mdl, fs, staged=False)
SCENARIOS COMPLETE: 66%|██████▌ | 78/118 [04:54<02:31, 3.78s/it]
Error at t=11.2 in scenario SingleFaultScenario(sequence={2.0: Injection(faults={'affect_dof': ['rf_propstuck']}, disturbances={})}, times=(2,), function='affect_dof', fault='rf_propstuck', rate=1.0416666666666668e-07, name='affect_dof_rf_propstuck_t2', time=2, phase='move')
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[22], line 1
----> 1 endresults, hists = propagate.fault_sample(mdl, fs, staged=False)
File ~\Documents\GitHub\fmdtools\fmdtools\sim\propagate.py:631, in fault_sample(mdl, fs, include_nominal, get_phasemap, **kwargs)
628 nomresult, nomhist, nomscen, c_mdl, t_end_nom = n_outs
629 scenlist = fs.scenarios()
--> 631 results, mdlhists = scenlist_helper(mdl,
632 scenlist,
633 c_mdl,
634 **kwargs,
635 nomhist=nomhist,
636 nomresult=nomresult)
638 if include_nominal:
639 process_nominal(mdlhists, nomhist, results, nomresult, t_end_nom, **kwargs)
File ~\Documents\GitHub\fmdtools\fmdtools\sim\propagate.py:803, in scenlist_helper(mdl, scenlist, c_mdl, **kwargs)
801 else:
802 mdl_i = c_mdl[0].new()
--> 803 ec, mh, t_end = exec_scen(mdl_i, scen, indiv_id=str(i), **kwargs)
804 results[name], mdlhists[name] = ec, mh
805 return results, mdlhists
File ~\Documents\GitHub\fmdtools\fmdtools\sim\propagate.py:839, in exec_scen(mdl, scen, save_args, indiv_id, **kwargs)
821 def exec_scen(mdl, scen, save_args={}, indiv_id='', **kwargs):
822 """
823 Executes a scenario and generates results and classifications given a model and
824 nominal model history.
(...)
837 :data:`sim_kwargs` for :func:`prop_one_scen`
838 """
--> 839 result, mdlhist, _, t_end, = prop_one_scen(mdl, scen, **kwargs)
840 save_helper(save_args, result, mdlhist, indiv_id=indiv_id, result_id=str(scen.name))
841 return result, mdlhist, t_end
File ~\Documents\GitHub\fmdtools\fmdtools\sim\propagate.py:1105, in prop_one_scen(mdl, scen, ctimes, nomhist, nomresult, **kwargs)
1103 disturbances = {}
1104 try:
-> 1105 mdl.propagate(t, fxnfaults, disturbances, run_stochastic=run_stochastic)
1106 except Exception as e:
1107 raise Exception("Error in scenario " + str(scen)) from e
File ~\Documents\GitHub\fmdtools\fmdtools\define\architecture\function.py:780, in FunctionArchitecture.propagate(self, time, fxnfaults, disturbances, run_stochastic)
778 # Step 2: Run Static Propagation Methods
779 try:
--> 780 self.prop_static(time, run_stochastic=run_stochastic)
781 except Exception as e:
782 raise Exception("Error in static propagation at time t=" + str(time)) from e
File ~\Documents\GitHub\fmdtools\fmdtools\define\architecture\function.py:809, in FunctionArchitecture.prop_static(self, time, run_stochastic)
806 flows_to_check = {*self.staticflows}
807 for fxnname in list(activefxns).copy():
808 # Update functions with new values, check to see if new faults or states
--> 809 oldmutables = self.fxns[fxnname].return_mutables()
810 self.fxns[fxnname]('static', time=time, run_stochastic=run_stochastic)
811 if oldmutables != self.fxns[fxnname].return_mutables():
File ~\Documents\GitHub\fmdtools\fmdtools\define\object\base.py:637, in BaseObject.return_mutables(self)
625 def return_mutables(self):
626 """
627 Return all mutable values in the block.
628
(...)
634 tuple of all mutable roles for the object.
635 """
636 return tuple([mut.return_mutables() if hasattr(mut, 'return_mutables')
--> 637 else mut for mut in self.find_mutables()])
File ~\Documents\GitHub\fmdtools\fmdtools\define\object\base.py:622, in BaseObject.find_mutables(self)
620 def find_mutables(self):
621 """Return list of mutable roles."""
--> 622 return [v for v in self.get_roles_as_dict(with_immutable=False).values()
623 if not ismethod(v)]
File ~\Documents\GitHub\fmdtools\fmdtools\define\object\base.py:501, in BaseObject.get_roles_as_dict(self, with_immutable, with_prefix, flex_prefixes, with_flex, no_flows, *roletypes, **kwargs)
498 if not non_flex_roletypes:
499 non_flex_roletypes = 'none'
--> 501 roles = self.get_roles(*non_flex_roletypes,
502 with_immutable=with_immutable,
503 no_flows=no_flows)
504 non_flex_roles = {role: getattr(self, role) for role in roles}
505 all_roles = {**flex_roles, **non_flex_roles}
File ~\Documents\GitHub\fmdtools\fmdtools\define\object\base.py:461, in BaseObject.get_roles(self, with_immutable, no_flows, *roletypes, **kwargs)
459 """Get all roles."""
460 roletypes = self.get_default_roletypes(*roletypes, no_flows=no_flows)
--> 461 return [role for roletype in roletypes
462 for role in getattr(self, roletype+'s', [])
463 if with_immutable or role not in self.immutable_roles]
File ~\Documents\GitHub\fmdtools\fmdtools\define\object\base.py:461, in <listcomp>(.0)
459 """Get all roles."""
460 roletypes = self.get_default_roletypes(*roletypes, no_flows=no_flows)
--> 461 return [role for roletype in roletypes
462 for role in getattr(self, roletype+'s', [])
463 if with_immutable or role not in self.immutable_roles]
KeyboardInterrupt:
[23]:
plot_env_with_traj_z(hists , mdl)
[24]:
statsfmea = tabulate.FMEA(endresults, fs,
weight_metrics=['rate'],
avg_metrics=['unsafe_flight_time', 'cost', 'repcost',
'landcost', 'body_strikes',
'head_strikes', 'property_restrictions'])
fmeatab = statsfmea.as_table(sort_by="cost")
fmeatab
[ ]: