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.

[15]:

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.

[16]:

from drone_mdl_urban import Drone

[17]:

mdl = Drone()

This is the model structure:

[18]:

mg = mdl.as_modelgraph()
fig, ax = mg.draw()
../../_images/examples_multirotor_Urban_Drone_Demo_6_0.png

We can also view the grid environment using its show methods:

[19]:

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)
../../_images/examples_multirotor_Urban_Drone_Demo_8_0.png

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:

[20]:

fig, ax = mdl.flows['environment'].c.show_z("height", collections=collections)
../../_images/examples_multirotor_Urban_Drone_Demo_10_0.png

Nominal Simulation

Below we show how this drone performs in the nominal scenario.

[21]:

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')
../../_images/examples_multirotor_Urban_Drone_Demo_12_0.png

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:

[22]:

fig, ax = hist_nom.plot_trajectories('dofs.s.x', 'dofs.s.y', 'dofs.s.z')
../../_images/examples_multirotor_Urban_Drone_Demo_14_0.png

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.

[23]:

from drone_mdl_urban import plot_env_with_traj, plot_env_with_traj_z

[24]:

plot_env_with_traj_z(hist_nom, mdl)

[24]:

(<Figure size 400x400 with 1 Axes>,
 <Axes3D: title={'center': 'trajectory'}, xlabel='x', ylabel='y', zlabel='z'>)
../../_images/examples_multirotor_Urban_Drone_Demo_17_1.png 
[25]:

plot_env_with_traj(hist_nom, mdl)

[25]:

(<Figure size 500x500 with 2 Axes>,
 <Axes: title={'center': 'trajectory'}, xlabel='x', ylabel='y'>)
../../_images/examples_multirotor_Urban_Drone_Demo_18_1.png

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:

[26]:

results_nom

[26]:

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:

[27]:

mdl.fxns['affect_dof'].get_faults(only_present=False)

[27]:

{'drone.fxns.affect_dof': {},
 'drone.fxns.affect_dof.ca': {},
 'drone.fxns.affect_dof.ca.comps.lf': {'ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases=(), disturbances=(), units='sim'),
  'ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases=(), disturbances=(), units='sim'),
  'openc': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim'),
  'propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'short': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim')},
 'drone.fxns.affect_dof.ca.comps.lr': {'ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases=(), disturbances=(), units='sim'),
  'ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases=(), disturbances=(), units='sim'),
  'openc': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim'),
  'propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'short': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim')},
 'drone.fxns.affect_dof.ca.comps.rf': {'ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases=(), disturbances=(), units='sim'),
  'ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases=(), disturbances=(), units='sim'),
  'openc': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim'),
  'propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'short': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim')},
 'drone.fxns.affect_dof.ca.comps.rr': {'ctlbreak': Fault(prob=2.0000000000000003e-06, cost=1000, phases=(), disturbances=(), units='sim'),
  'ctldn': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'ctlup': Fault(prob=2.0000000000000003e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechbreak': Fault(prob=1.0000000000000002e-06, cost=500, phases=(), disturbances=(), units='sim'),
  'mechfriction': Fault(prob=5.000000000000001e-07, cost=500, phases=(), disturbances=(), units='sim'),
  'openc': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim'),
  'propbreak': Fault(prob=3.0000000000000004e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propstuck': Fault(prob=2.0000000000000002e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'propwarp': Fault(prob=1.0000000000000001e-07, cost=200, phases=(), disturbances=(), units='sim'),
  'short': Fault(prob=1.0000000000000002e-06, cost=200, phases=(), disturbances=(), units='sim')}}

[30]:

results_fault, hist_fault =propagate.one_fault(mdl, "affect_dof.ca.comps.lr", "mechbreak", time=2.0)

[31]:

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')
../../_images/examples_multirotor_Urban_Drone_Demo_24_0.png 
[32]:

plot_env_with_traj_z(hist_fault, mdl)

[32]:

(<Figure size 400x400 with 1 Axes>,
 <Axes3D: title={'center': 'trajectory'}, xlabel='x', ylabel='y', zlabel='z'>)
../../_images/examples_multirotor_Urban_Drone_Demo_25_1.png 
[33]:

plot_env_with_traj(hist_fault, mdl)

[33]:

(<Figure size 500x500 with 2 Axes>,
 <Axes: title={'center': 'trajectory'}, xlabel='x', ylabel='y'>)
../../_images/examples_multirotor_Urban_Drone_Demo_26_1.png 
[34]:

results_fault

[34]:

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:

[35]:

from fmdtools.analyze.phases import PhaseMap, from_hist
phasemaps = from_hist(hist_nom, fxn_modephases=[])
phasemaps

[35]:

{'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]}, {})}

[36]:

from drone_mdl_urban import make_move_quad
move_quad=make_move_quad(hist_nom, phasemaps['plan_path'].phases['move'])
move_quad

[36]:

{'samp': 'quadrature',
 'quad': {'nodes': [-0.20000000000000018, 1.0], 'weights': [1.0, 0.92]}}

[37]:

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

[37]:

FaultSample of scenarios:
 - drone_fxns_manage_health_lostfunction_t2
 - drone_fxns_manage_health_lostfunction_t2
 - drone_fxns_store_ee_lowcharge_t2
 - drone_fxns_store_ee_lowcharge_t2
 - drone_fxns_store_ee_nocharge_t2
 - drone_fxns_store_ee_nocharge_t2
 - drone_fxns_store_ee_ca_comps_s1p1_break_t2
 - drone_fxns_store_ee_ca_comps_s1p1_break_t2
 - drone_fxns_store_ee_ca_comps_s1p1_degr_t2
 - drone_fxns_store_ee_ca_comps_s1p1_degr_t2
 - ... (118 total)

[38]:

endresults, hists = propagate.fault_sample(mdl, fs, staged=False)

SCENARIOS COMPLETE: 100%|██████████| 118/118 [01:32<00:00,  1.27it/s]

[39]:

plot_env_with_traj_z(hists , mdl)

[39]:

(<Figure size 400x400 with 1 Axes>,
 <Axes3D: title={'center': 'trajectory'}, xlabel='x', ylabel='y', zlabel='z'>)
../../_images/examples_multirotor_Urban_Drone_Demo_33_1.png 
[ ]:

statsfmea = tabulate.FMEA(endresults, fs,
                          average_metric=['rate', 'unsafe_flight_time', 'cost', 'repcost',
                                            'landcost', 'body_strikes',
                                            'head_strikes', 'property_restrictions'],
                          rates='rate')
fmeatab = statsfmea.as_table(sort_by="cost")
fmeatab
average_rate average_unsafe_flight_time average_cost average_repcost average_landcost average_body_strikes average_head_strikes average_property_restrictions
drone.fxns.plan_path.ca.comps.vision undesired_detection 5.989583e-03 0.0 0.0 0.0 0.0 0.0 0.0 0.0
lack_of_detection 5.989583e-03 281.0 28100.0 0.0 0.0 0.0 0.0 0.0
drone.fxns.store_ee lowcharge 6.708333e-06 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
drone.fxns.dist_ee degr 2.395833e-06 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.store_ee nocharge 1.916667e-06 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
drone.fxns.dist_ee short 1.437500e-06 281.0 29600.0 1500.0 0.0 0.0 0.0 0.0
drone.fxns.plan_path degloc 1.150000e-06 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.ctl_dof degctl 1.150000e-06 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr ctlup 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf ctldn 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr ctlbreak 9.583333e-07 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr ctlbreak 9.583333e-07 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf ctlbreak 9.583333e-07 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr ctldn 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr ctldn 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf ctlup 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf ctlup 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf ctlbreak 9.583333e-07 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.dist_ee break 9.583333e-07 281.0 29600.0 1500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr ctlup 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf ctldn 9.583333e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf mechbreak 4.791667e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr openc 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf short 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr mechbreak 4.791667e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr mechbreak 4.791667e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf mechbreak 4.791667e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
openc 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
short 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr openc 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr short 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf openc 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr short 4.791667e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.ctl_dof noctl 2.875000e-07 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.plan_path noloc 2.875000e-07 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf mechfriction 2.395833e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr mechfriction 2.395833e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr mechfriction 2.395833e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf mechfriction 2.395833e-07 281.0 28600.0 500.0 0.0 0.0 0.0 0.0
drone.fxns.store_ee.ca.comps.s1p1 nocharge 1.437500e-07 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf propbreak 1.437500e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr propbreak 1.437500e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr propbreak 1.437500e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf propbreak 1.437500e-07 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.store_ee.ca.comps.s1p1 lowcharge 9.583333e-08 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr propstuck 9.583333e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf propstuck 9.583333e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf propstuck 9.583333e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr propstuck 9.583333e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.store_ee.ca.comps.s1p1 short 7.187500e-08 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
break 7.187500e-08 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
degr 7.187500e-08 281.0 28200.0 100.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rr propwarp 4.791667e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.rf propwarp 4.791667e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lr propwarp 4.791667e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.affect_dof.ca.comps.lf propwarp 4.791667e-08 281.0 28300.0 200.0 0.0 0.0 0.0 0.0
drone.fxns.manage_health lostfunction 7.187500e-09 281.0 29100.0 1000.0 0.0 0.0 0.0 0.0
drone.fxns.hold_payload deform 2.875000e-09 281.0 29600.0 1500.0 0.0 0.0 0.0 0.0
break 7.187500e-10 281.0 29600.0 1500.0 0.0 0.0 0.0 0.0 
[ ]: