Updated demo, made pyscurve optional in setup.cfg (09aaadc8) · Commits · Ibrahim Ahmed / multirotor

Demo.ipynb

+88 −122

Changes for Demo.ipynb: 88 added lines, 122 removed lines.
Original line number	Diff line number	Diff line
		%% Cell type:markdown id:d001f534 tags:

		#### Notebook setup

		%% Cell type:code id:be76bd53 tags:

		``` python
		%pip install -e .
		```

		%% Cell type:code id:5a776232 tags:

		``` python
		%reload_ext autoreload
		%autoreload 2
		%matplotlib inline
		import time
		import warnings
		import os, sys
		from copy import deepcopy
		from types import SimpleNamespace

		import matplotlib.pyplot as plt
		import gym
		import numpy as np
		from tqdm.auto import tqdm, trange

		from multirotor.helpers import control_allocation_matrix, DataLog
		from multirotor.vehicle import MotorParams, VehicleParams, PropellerParams, SimulationParams, BatteryParams
		from multirotor.controller import (
		PosController, VelController,
		AttController, RateController,
		AltController, AltRateController,
		Controller
		)
		from multirotor.simulation import Multirotor, Propeller, Motor, Battery
		from multirotor.coords import body_to_inertial, inertial_to_body, direction_cosine_matrix, euler_to_angular_rate
		from multirotor.env import SpeedsMultirotorEnv as LocalOctorotor
		from multirotor.trajectories import Trajectory, GuidedTrajectory
		from multirotor.env import SpeedsMultirotorEnv, DynamicsMultirotorEnv
		from multirotor.trajectories import Trajectory
		from multirotor.visualize import plot_datalog
		```

		%% Cell type:code id:421b284b tags:

		``` python
		# Plotting/display parameters
		# https://stackoverflow.com/a/21009774/4591810
		float_formatter = "{:.3f}".format
		np.set_printoptions(formatter={'float_kind':float_formatter})

		SMALL_SIZE = 16
		MEDIUM_SIZE = 16
		BIGGER_SIZE = 20

		plt.rc('font', size=SMALL_SIZE) # controls default text sizes
		plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
		plt.rc('axes', labelsize=BIGGER_SIZE, titlesize=BIGGER_SIZE) # fontsize of the x and y labels
		plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
		plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
		plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
		plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
		```

		%% Cell type:markdown id:9a31149f tags:

		### Parameters

		%% Cell type:code id:3c2c1785 tags:

		``` python
		22.2 / 0.0347
		```

		%% Cell type:code id:1ba2a2fe tags:

		``` python
		# Tarot T18 params
		bp = BatteryParams(max_voltage=22.2)
		mp = MotorParams(
		moment_of_inertia=5e-5,
		# resistance=0.27,
		resistance=0.081,
		k_emf=0.0265,
		# k_motor=0.0932,
		speed_voltage_scaling=0.0347,
		speed_voltage_scaling= 0.0347,
		max_current=38.
		)
		pp = PropellerParams(
		moment_of_inertia=1.86e-6,
		use_thrust_constant=True,
		k_thrust=9.8419e-05, # 18-inch propeller
		# k_thrust=5.28847e-05, # 15 inch propeller
		k_drag=1.8503e-06, # 18-inch propeller
		# k_drag=1.34545e-06, # 15-inch propeller
		motor=mp
		)
		vp = VehicleParams(
		propellers=[pp] * 8,
		battery=bp,
		# angles in 45 deg increments, rotated to align with
		# model setup in gazebo sim (not part of this repo)
		angles=np.linspace(0, -2np.pi, num=8, endpoint=False) + 0.375 np.pi,
		distances=np.ones(8) * 0.635,
		clockwise=[-1,1,-1,1,-1,1,-1,1],
		mass=10.66,
		inertia_matrix=np.asarray([
		[0.2206, 0, 0],
		[0, 0.2206, 0.],
		[0, 0, 0.4238]
		])
		)
		sp = SimulationParams(dt=0.01, g=9.81)
		```

		%% Cell type:markdown id:17f3b34d tags:

		### Multirotor

		%% Cell type:markdown id:0505e4f6 tags:

		#### Motor

		%% Cell type:code id:b3128c10 tags:

		``` python
		%matplotlib inline
		# Plot motor speeds as a function of time and input voltage signal
		plt.figure(figsize=(8,8))
		motor = Motor(mp, sp)
		for signal in [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]:
		for vsignal in [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]:
		speeds = []
		motor.reset()
		speed = vsignal / mp.speed_voltage_scaling
		for i in range(200):
		speeds.append(motor.step(signal))
		plt.plot(speeds, label='%dV' % signal)
		speeds.append(motor.step(speed))
		plt.plot(speeds, label='%d rad/s' % speed)
		plt.legend(ncol=2)
		plt.ylabel('Speed rad/s')
		plt.xlabel('Time /ms')
		```

		%% Cell type:code id:42ebfb97 tags:

		``` python
		from multirotor.helpers import learn_speed_voltage_scaling

		def make_motor_fn(params, sp):
		from copy import deepcopy
		params = deepcopy(params)
		params.speed_voltage_scaling = 1.
		def motor_step(signal):
		m = Motor(params, sp)
		for i in range(100):
		s = m.step(signal)
		return s
		return motor_step

		print('Voltage = %.5f * speed' % (learn_speed_voltage_scaling(make_motor_fn(mp, sp))))
		```

		%% Cell type:markdown id:6fb2f5c9 tags:

		#### Propeller

		%% Cell type:code id:1669687b tags:

		``` python
		%matplotlib inline
		# Plot propeller speed by numerically solving the thrust equation,
		# if accurate propeller measurements are given in params
		pp_ = deepcopy(pp)
		pp_.use_thrust_constant = False # Set to true to just use k_thrust
		prop = Propeller(pp_, sp)
		plt.figure(figsize=(8,8))
		speeds = np.linspace(0, 600, num=100)
		for a in np.linspace(0, 10, 10, endpoint=False):
		thrusts = []
		for s in speeds:
		thrusts.append(prop.thrust(s, np.asarray([0, 0, a])))
		plt.plot(speeds, thrusts, label='%.1f m/s' % a)
		plt.xlabel('Speed rad/s')
		plt.ylabel('Thrust /N')
		plt.title('Thrust with airspeed')
		plt.legend(ncol=2)
		```

		%% Cell type:markdown id:e7bb489d tags:

		#### Vehicle

		%% Cell type:code id:ecc23e14 tags:

		``` python
		# Combine propeller/motor/vehicle to get vehicle.
		# Take off simulation
		m = Multirotor(vp, sp)
		log = DataLog(vehicle=m) # convenient logging class
		m.reset()
		m.state *= 0 # set to zero, reset() sets random values
		action = m.allocate_control(
		thrust=m.weight * 1.1,
		torques=np.asarray([0, 0, 0])
		)
		for i in range(500):
		m.step_speeds(action)
		log.log()
		log.done_logging()
		plt.plot(log.z)
		```

		%% Cell type:markdown id:2c6535f4 tags:

		### PID Controller

		%% Cell type:code id:33984887 tags:

		``` python
		# From PID parameters file
		def get_controller(m: Multirotor, max_velocity=5., max_acceleration=3.):
		assert m.simulation.dt <= 0.1, 'Simulation time step too large.'
		pos = PosController(
		1.0, 0., 0., 1., vehicle=m,
		max_velocity=max_velocity, max_acceleration=max_acceleration
		max_velocity=max_velocity, max_acceleration=max_acceleration,
		square_root_scaling=True, leashing=True
		)
		vel = VelController(
		2.0, 1.0, 0.5, 1000., vehicle=m)
		2.0, 1.0, 0.5,
		max_err_i=max_acceleration,
		max_tilt=np.pi/12,
		vehicle=m)
		att = AttController(
		[2.6875, 4.5, 4.5],
		0, 0.,
		1., vehicle=m)
		max_err_i=1.,
		vehicle=m)
		rat = RateController(
		[4., 4., 4.],
		0, 0, # purely P control
		# [0.1655, 0.1655, 0.5],
		# [0.135, 0.135, 0.018],
		# [0.01234, 0.01234, 0.],
		[0.5,0.5,0.5],
		0, 0,
		max_err_i=0.5,
		max_acceleration=1.,
		vehicle=m)
		alt = AltController(
		1, 0, 0,
		1, vehicle=m)
		max_err_i=1, vehicle=m,
		max_velocity=max_velocity)
		alt_rate = AltRateController(
		5, 0, 0,
		1, vehicle=m)
		max_err_i=1, vehicle=m)
		ctrl = Controller(
		pos, vel, att, rat, alt, alt_rate,
		interval_p=0.1, interval_a=0.01, interval_z=0.1
		period_p=0.1, period_a=0.01, period_z=0.1
		)
		return ctrl
		```

		%% Cell type:code id:d21d73fe tags:

		``` python
		%matplotlib inline
		m = Multirotor(vp, sp)
		ctrl = get_controller(m)
		log = DataLog(vehicle=m, controller=ctrl)
		for i in range(500):
		action = ctrl.step((0.01,0.1,1,0))
		# no allocation or motor simulation, for which we first need to
		# m.step_speeds(m.allocate_control(action[0], action[3:])
		# Instead, requested dynamics are fulfilled:
		dynamics = np.zeros(6, m.dtype)
		dynamics[2] = action[0]
		dynamics[3:] = action[1:]
		# no allocation or motor simulation,
		# requested dynamics are fulfilled:
		m.step_dynamics(dynamics)
		log.log()
		log.done_logging()

		plt.plot(log.actions[:,0], ls=':', label='thrust')
		lines = plt.gca().lines
		plt.twinx()
		for s, axis in zip(log.actions.T[1:], ('x','y','z')):
		plt.plot(s, label=axis + '-torque')
		plt.legend(handles=plt.gca().lines + lines)
		```

		%% Cell type:markdown id:cfbc9c25 tags:

		#### Attitude Angle Controller

		%% Cell type:code id:f012e1f8 tags:

		``` python
		m = Multirotor(vp, sp)
		fz = m.weight
		att = get_controller(m).ctrl_a
		log = DataLog(vehicle=m, controller=att, other_vars=('err',))
		for i in range(5000):
		ref = np.asarray([np.pi/18, 0, 0])
		# action is prescribed euler rate
		action = att.step(ref, m.orientation)
		# action = np.clip(action, a_min=-0.1, a_max=0.1)
		m.step_dynamics(np.asarray([0, 0, 0, *action]))
		log.log(err=att.err_p[0])
		log._actions[-1] = action
		log.done_logging()

		plt.plot(log.roll * 180 / np.pi)
		plt.twinx()
		plt.plot(log.actions[:,0], ls=':', label='Rate rad/s')
		```

		%% Cell type:markdown id:fdbd88c1 tags:

		#### Attitude Rate Controller

		%% Cell type:code id:9cc3b317 tags:

		``` python
		m = Multirotor(vp, sp)
		fz = m.weight
		ctrl = get_controller(m)
		rat = ctrl.ctrl_r
		att = ctrl.ctrl_a
		log = DataLog(vehicle=m, controller=rat, other_vars=('err',))
		for i in range(200):
		ref = np.asarray([np.pi/18, np.pi/12, 0])
		rate = att.step(ref, m.orientation, m.simulation.dt)
		torque = rat.step(rate, m.euler_rate, m.simulation.dt)
		action = np.clip(torque, a_min=-0.1, a_max=0.1)
		m.step_dynamics(np.asarray([0, 0, 0, *action]))
		log.log(err=rat.err_p[0])
		log._actions[-1] = action
		log.done_logging()

		plt.plot(log.roll * 180 / np.pi, c='r', label='roll')
		plt.plot(log.pitch * 180 / np.pi, c='g', label='pitch')
		plt.plot(log.yaw * 180 / np.pi, c='b', label='yaw')
		plt.ylabel('Orientation /deg')
		plt.legend()
		plt.twinx()
		plt.plot(log.actions[:,0], ls=':', c='r')
		plt.plot(log.actions[:,1], ls=':', c='g')
		plt.plot(log.actions[:,2], ls=':', c='b')
		plt.ylabel('Torque / Nm')
		plt.title('Ref orientation' + str(ref))
		```

		%% Cell type:code id:a7ad3618 tags:

		``` python
		m.params.inertia_matrix
		```

		%% Cell type:code id:f996d427 tags:

		``` python
		att.err
		```

		%% Cell type:code id:8b3b2181 tags:

		``` python
		rat.err
		```

		%% Cell type:code id:d5b77d12 tags:

		``` python
		action
		```

		%% Cell type:code id:3b9a3eba tags:

		``` python
		euler_to_angular_rate(
		m.euler_rate,
		m.orientation
		)
		```

		%% Cell type:markdown id:6370aa80 tags:

		#### Altitude Controller

		%% Cell type:code id:5fcd0d24 tags:

		``` python
		m = Multirotor(vp, sp)
		ctrl = get_controller(m)
		alt = ctrl.ctrl_z
		alt_rate = ctrl.ctrl_vz
		log = DataLog(vehicle=m, controller=alt, other_vars=('thrust',))
		log = DataLog(vehicle=m, other_vars=('thrust',))
		for i in range(5000):
		ref = np.asarray([1.])
		rate = alt.step(ref, m.position[2:])
		action = alt_rate.step(rate, m.world_velocity[2:])
		rate = alt.step(ref, m.position[2:], dt=0.1)
		action = alt_rate.step(rate, m.inertial_velocity[2:], dt=0.1)
		action = np.clip(action, a_min=-2m.weight, a_max=2m.weight)
		m.step_dynamics(np.asarray([0, 0, action[0], 0,0,0]))
		log.log(thrust=action)
		#log._actions[-1] = action
		log.done_logging()

		plt.plot(log.actions.squeeze())
		plt.plot(log.thrust.squeeze())
		plt.twinx()
		plt.plot(log.z, ls=':')
		```

		%% Cell type:markdown id:f4278c17 tags:

		#### Position Controller

		%% Cell type:code id:91ae8919 tags:

		``` python
		m = Multirotor(vp, sp)
		ctrl = get_controller(m)
		pos = ctrl.ctrl_p
		vel = ctrl.ctrl_v
		rat = ctrl.ctrl_r
		att = ctrl.ctrl_a
		log = DataLog(vehicle=m, controller=pos, other_vars=('err', 'att_actions'))
		log = DataLog(vehicle=m, other_vars=('err', 'att_actions'))
		for i in range(5000):
		ref = np.asarray([1.,0.])
		velocity = pos.step(ref, m.position[:2])
		angles = vel.step(velocity, m.velocity[:2])[::-1]
		rate = att.step(np.asarray([*angles, 0]), m.orientation)
		action = rat.step(rate, m.euler_rate)

		velocity = pos.step(ref, m.position[:2], dt=0.1)
		angles = vel.step(velocity, m.velocity[:2], dt=0.1)[::-1]
		rate = att.step(np.asarray([*angles, 0]), m.orientation, dt=0.01)
		action = rat.step(rate, m.euler_rate, dt=0.01)

		action = np.clip(action, a_min=-0.1, a_max=0.1)
		m.step_dynamics(np.asarray([0, 0, m.weight, *action]))
		log.log(err=pos.err_p[0], att_actions=action)
		log.done_logging()

		plt.plot(log.position[:,0])
		plt.plot(log.err)
		# plt.plot(log.position[:,1])
		plt.twinx()
		plt.plot(log.actions[:,0] * 180 / np.pi, ls=':')
		plt.plot(log.att_actions[:,0] * 180 / np.pi, ls=':')
		plt.plot(log.pitch * 180 / np.pi, ls='-.')
		# plt.plot(log.actions[:,0] * 180 / np.pi, ls=':')
		```

		%% Cell type:markdown id:565cf2c9 tags:

		#### Parameter search #TODO

		%% Cell type:code id:11013f70 tags:

		``` python
		# TODO
		```

		%% Cell type:markdown id:9eb34672 tags:

		### Simulation

		%% Cell type:code id:b60ef8d3 tags:

		``` python
		# NASA flight test
		# wp = np.asarray([
		# [0.0, 0.0, 30.0],
		# [164.0146725649829, -0.019177722744643688, 30.0],
		# [165.6418055187678, 111.5351051245816, 30.0],
		# [127.3337449710234, 165.73576059611514, 30.0],
		# [-187.28170707810204, 170.33217775914818, 45.0],
		# [-192.03130502498243, 106.30660058604553, 45.0],
		# [115.89920266153058, 100.8644210617058, 30.0],
		# [114.81859536317643, 26.80923518165946, 30.0],
		# [-21.459931490011513, 32.60508110653609, 30.0]
		# ])
		# wp = np.asarray([
		# [0.0, 0.0, 100.0],
		# [300, 0.0, 100.0],
		# [300, 200, 100.0],
		# [250, 300, 100.0],
		# [-350, 340, 120.0],
		# [-350, 250, 120.0],
		# [300, 250, 100.0],
		# [300, 150, 100.0],
		# [-50, 100, 100.0]
		# ])

		wp = np.asarray([
		[0.0, 0.0, 30.0],
		[164.0146725649829, -0.019177722744643688, 30.0],
		[165.6418055187678, 111.5351051245816, 30.0],
		[127.3337449710234, 165.73576059611514, 30.0],
		[-187.28170707810204, 170.33217775914818, 45.0],
		[-192.03130502498243, 106.30660058604553, 45.0],
		[115.89920266153058, 100.8644210617058, 30.0],
		[114.81859536317643, 26.80923518165946, 30.0],
		[-21.459931490011513, 32.60508110653609, 30.0]
		[0,10,0],
		[10,10,0],
		[10,0,0],
		[0,0,0]
		])
		```

		%% Cell type:code id:f5f52df4 tags:

		``` python
		def wind(t, m, nominal=False):
		if nominal:
		return np.zeros(3)
		w_inertial = np.asarray([5 * np.sin(t * 2 * np.pi / 4000), 0, 0])
		dcm = direction_cosine_matrix(*m.orientation)
		return inertial_to_body(w_inertial, dcm)
		```

		%% Cell type:code id:a98de44f tags:

		``` python
		def run_sim(env, traj, steps=60_000, disturbance=None):
		ctrl = get_controller(env.vehicle, max_velocity=5.)
		ctrl = get_controller(env.vehicle, max_velocity=3.)

		log = DataLog(env.vehicle, ctrl,
		other_vars=('speeds','target', 'alloc_errs', 'att_err',
		'rate_target', 'att_target',
		'leash', 'currents', 'voltages'))
		other_vars=('currents', 'voltages'))
		disturb_force, disturb_torque = 0., 0
		for i, (pos, feed_forward_vel) in tqdm(
		enumerate(traj), leave=False, total=steps
		):
		if i==steps: break
		# Generate reference for controller
		ref = np.asarray([*pos, 0.])
		# Get prescribed dynamics for system as thrust and torques
		dynamics = ctrl.step(ref, feed_forward_velocity=feed_forward_vel)
		thrust, torques = dynamics[0], dynamics[1:]
		# Allocate control: Convert dynamics into motor rad/s
		action = env.vehicle.allocate_control(thrust, torques)
		# get any disturbances
		if disturbance is not None:
		disturb_force, disturb_torque = disturbance(i, env.vehicle)
		# Send speeds to environment
		state, *_ = env.step(
		action, disturb_forces=disturb_force, disturb_torques=disturb_torque
		)
		alloc_errs = np.asarray([thrust, torques]) - env.vehicle.alloc @ action*2

		log.log(speeds=action, target=pos, alloc_errs=alloc_errs,
		leash=ctrl.ctrl_p.leash,
		att_err=ctrl.ctrl_a.err,
		att_target = ctrl.ctrl_v.action[::-1],
		rate_target=ctrl.ctrl_a.action,
		currents=[p.motor.current for p in env.vehicle.propellers],
		log.log(currents=[p.motor.current for p in env.vehicle.propellers],
		voltages=[p.motor.voltage for p in env.vehicle.propellers])

		if np.any(np.abs(env.vehicle.orientation[:2]) > np.pi/6): break

		log.done_logging()
		return log
		```

		%% Cell type:code id:62689b31 tags:

		``` python
		env = LocalOctorotor(vehicle=Multirotor(vp, sp))
		env = SpeedsMultirotorEnv(vehicle=Multirotor(vp, sp)) # step() takes speeds

		# waypoints = [[0,50,2], [50,50,2], [50,0,2], [0,0,2]]
		waypoints = wp
		# traj = GuidedTrajectory(env.vehicle, waypoints, proximity=2)
		traj = Trajectory(env.vehicle, waypoints, proximity=2, resolution=None)
		log = run_sim(env, traj, steps=60_000, disturbance=None)
		```

		%% Cell type:code id:2d1b0058 tags:

		``` python
		# Currents
		plt.plot(log.currents, ls=':')
		plt.ylabel('Motor current /A')
		plt.xlabel('Time /ms')
		plt.title('Individual motor currents')
		```

		%% Cell type:code id:f0aa49e1 tags:

		``` python
		# Voltages
		plt.plot(log.voltages, ls=':')
		plt.ylim(0, 30)
		plt.ylabel('Motor voltage /A')
		plt.xlabel('Time /ms')
		plt.title('Voltages')
		```

		%% Cell type:code id:81a8043e tags:
		%% Cell type:code id:665f4799 tags:

		``` python
		%matplotlib inline
		plt.figure(figsize=(21,10.5))
		plot_grid = (3,3)
		plt.subplot(*plot_grid,1)

		n = len(log)

		plt.plot(log.x, label='x', c='r')
		plt.plot(log.target[:, 0], c='r', ls=':')
		plt.plot(log.y, label='y', c='g')
		plt.plot(log.target[:, 1], c='g', ls=':')
		plt.plot(log.z, label='z', c='b')
		lines = plt.gca().lines[::2]
		plt.ylabel('Position /m')
		plt.twinx()
		plt.plot(log.roll * (180 / np.pi), label='roll', c='c', ls=':')
		plt.plot(log.pitch * (180 / np.pi), label='pitch', c='m', ls=':')
		plt.plot(log.yaw * (180 / np.pi), label='yaw', c='y', ls=':')
		plt.ylabel('Orientation /deg')
		plt.legend(handles=plt.gca().lines + lines, ncol=2)
		plt.title('Position and Orientation')

		plt.subplot(*plot_grid,2)
		for i in range(log.speeds.shape[1]):
		l, = plt.plot(log.speeds[:,i], label='prop %d' % i)
		# plt.plot(speeds[:,i], c=l.get_c())
		lines = plt.gca().lines
		plt.legend(handles=lines, ncol=2)
		plt.title('Motor speeds /RPM')


		plt.subplot(*plot_grid,3)
		v_world = np.zeros_like(log.velocity)
		for i, (v, o) in enumerate(zip(log.velocity, log.orientation)):
		dcm = direction_cosine_matrix(*o)
		v_world[i] = body_to_inertial(v, dcm)
		for i, c, a in zip(range(3), 'rgb', 'xyz'):
		plt.plot(v_world[:,i], label='Velocity %s' % a, c=c)
		# plt.plot(velocities[:,i], label='Velocity %s' % a, c=c)
		plt.legend()
		plt.title('Velocities')

		plt.subplot(*plot_grid,4)
		plt.title('Controller allocated dynamics')
		l = plt.plot(log.actions[:,0], label='Ctrl Thrust')
		plt.ylabel('Force /N')
		plt.twinx()
		for i, c, a in zip(range(3), 'rgb', 'xyz'):
		plt.plot(log.actions[:,1+i], label='Ctrl Torque %s' % a, c=c)
		plt.ylabel('Torque /Nm')
		plt.legend(handles=plt.gca().lines + l, ncol=2)

		plt.subplot(*plot_grid,5)
		lines = plt.plot(log.alloc_errs[:, 0], label='Thrust err', c='b')
		plt.ylabel('Thrust /N')
		plt.twinx()
		plt.plot(log.alloc_errs[:, 1], label='Torque x err', ls=':')
		plt.plot(log.alloc_errs[:, 2], label='Torque y err', ls=':')
		plt.plot(log.alloc_errs[:, 3], label='Torque z err', ls=':')
		plt.legend(handles = plt.gca().lines + lines, ncol=2)
		plt.ylabel('Torque /Nm')
		plt.title('Allocation Errors')

		plt.subplot(*plot_grid,6)
		plt.plot(log.target[:,0], log.target[:,1], label='Prescribed traj')
		plt.plot(log.x, log.y, label='Actual traj', ls=':')
		plt.gca().set_aspect('equal', 'box')
		plt.title('XY positions /m')
		plt.legend()

		plt.tight_layout()
		# PLot positions, velocities, prescribed dynamics
		plot_datalog(log)
		```

		%% Cell type:code id:d9d99572 tags:

		``` python
		# 3D plot of trajectory
		%matplotlib notebook
		fig = plt.figure()
		xlim = ylim = zlim = (np.min(log.position), np.max(log.position))
		ax = fig.add_subplot(projection='3d', xlim=xlim, ylim=ylim, zlim=zlim)
		ax.set_xlabel('x')
		ax.set_ylabel('y')
		ax.set_zlabel('z')
		ax.plot(log.x, log.y, log.z)
		```

multirotor/helpers.py

+7 −8

Changes for multirotor/helpers.py: 7 added lines, 8 removed lines.
Original line number	Diff line number	Diff line
		@@ -232,11 +232,8 @@ def control_allocation_matrix(params: VehicleParams) -> Tuple[np.ndarray, np.nda
		def get_vehicle_ability(
		vp: VehicleParams, sp: SimulationParams,
		max_tilt: float=np.pi/12,
		max_angular_acc: float=5,
		max_rads: float=600
		):
		alloc, alloc_inverse = control_allocation_matrix(vp)
		I_roll, I_pitch, I_yaw = vp.inertia_matrix.diagonal()
		n = len(vp.propellers)

		thrusts = [p.k_thrust * max_rads**2 for p in vp.propellers]
		@@ -316,11 +313,6 @@ class DataLog:
		controller : Controller
		The controller to track.
		"""
		self._states_names = ('x','y','z',
		'vx','vy','vz',
		'roll','pitch','yaw',
		'xrate', 'yrate', 'zrate')
		self._action_names = ('thrust', 'torque_x', 'torque_y', 'torque_z')
		self._arrayed = False
		self._states = []
		self.states = None
		@@ -329,6 +321,13 @@ class DataLog:
		self.times = None
		self._times = []
		self._args = () if other_vars is None else other_vars
		own_vars = ('arrayed', 'states', 'actions', 'times', 'args', 'target')
		not_allowed = [v for v in own_vars if v in self._args]
		if len(not_allowed) > 0:
		raise AttributeError(
		('The following `other_vars` are not allowed since they are attributes '
		', '.join(not_allowed))
		)
		for arg in self._args:
		setattr(self, arg, None)
		setattr(self, '_' + str(arg), [])

multirotor/simulation.py

+16 −13

Changes for multirotor/simulation.py: 16 added lines, 13 removed lines.
Original line number	Diff line number	Diff line
		@@ -189,7 +189,7 @@ class Motor:
		Parameters
		----------
		u : float
		Voltage signal.
		Speed signal (rad/s).
		max_voltage : float, optional
		The maximum voltage supply from power source. By default infinite.

		@@ -267,9 +267,7 @@ class Multirotor:
		self._dxdt = None
		self.dxdt_decimals = max(1, 1 - int(np.log10(self.simulation.dt)))

		self.propellers: List[Propeller] = None
		self.propeller_vectors: np.ndarray = None
		self.propellers = []
		self.propellers: List[Propeller] = []
		for params in self.params.propellers:
		self.propellers.append(Propeller(params, self.simulation))

		@@ -302,6 +300,7 @@ class Multirotor:
		self.alloc = self.alloc.astype(self.dtype)
		self.params.propeller_vectors = self.params.propeller_vectors.astype(self.dtype)
		self.alloc_inverse = self.alloc_inverse.astype(self.dtype)
		self.params.inertia_matrix = self.params.inertia_matrix.astype(self.dtype)
		self.params.inertia_matrix_inverse = self.params.inertia_matrix_inverse.astype(self.dtype)
		self.state = np.zeros(12, dtype=self.dtype)
		self._dxdt = np.zeros_like(self.state)
		@@ -450,10 +449,10 @@ class Multirotor:
		# Do not need to get forces/torques on body, since the action array
		# already is a 6d vector of forces/torques.
		# forces, torques = self.get_forces_torques(u, x)
		xdot = apply_forces_torques(
		u[:3], u[3:], x, self.simulation.g,
		dxdt = apply_forces_torques(
		u[:3], u[3:], x.astype(self.dtype), self.simulation.g,
		self.params.mass, self.params.inertia_matrix, self.params.inertia_matrix_inverse)
		return np.around(xdot, self.dxdt_decimals)
		return np.around(dxdt, self.dxdt_decimals)


		def dxdt_speeds(
		@@ -487,10 +486,12 @@ class Multirotor:
		# same state by the odeint() function, and the results should be consistent.
		forces, torques = self.get_forces_torques(
		u, x)
		xdot = apply_forces_torques(
		forces+disturb_forces, torques+disturb_torques, x, self.simulation.g,
		# print('dxdt-x', self.t // self.simulation.dt, x.dtype)
		dxdt = apply_forces_torques(
		forces+disturb_forces, torques+disturb_torques, x.astype(self.dtype), self.simulation.g,
		self.params.mass, self.params.inertia_matrix, self.params.inertia_matrix_inverse)
		return np.around(xdot, self.dxdt_decimals)
		# print('dxdt', self.t // self.simulation.dt, dxdt.dtype)
		return np.around(dxdt, self.dxdt_decimals)


		def step_dynamics(self, u: np.ndarray) -> np.ndarray:
		@@ -516,7 +517,7 @@ class Multirotor:
		args=(u,),
		rtol=1e-4, atol=1e-4, tfirst=True
		)[-1]
		self.state = np.around(self.state, 4)
		self.state = np.around(self.state, 4).astype(self.dtype)
		# TODO: inverse solve for speed = forces to set propeller speeds
		return self.state

		@@ -550,12 +551,14 @@ class Multirotor:
		t=self.t, x=self.state, u=u,
		disturb_forces=disturb_forces, disturb_torques=disturb_torques
		)
		# print('pre-x', self.t // self.simulation.dt, self.state.dtype)
		self.state = odeint(
		self.dxdt_speeds, self.state, (0, self.simulation.dt),
		args=(u, disturb_forces, disturb_torques),
		rtol=1e-4, atol=1e-4, tfirst=True
		)[-1]
		self.state = np.around(self.state, 4)
		self.state = np.around(self.state, 4).astype(self.dtype)
		# print('post-x', self.t // self.simulation.dt, self.state.dtype)
		for u_, prop in zip(u, self.propellers):
		prop.step(u_, max_voltage=self.battery.voltage)
		self.battery.step()
		@@ -580,7 +583,7 @@ class Multirotor:
		The prescribed propeller speeds (rad /s)
		"""
		# TODO: njit it? np.linalg.lstsq can be compiled
		vec = np.asarray([thrust, *torques])
		vec = np.asarray([thrust, *torques], self.dtype)
		# return np.sqrt(np.linalg.lstsq(self.alloc, vec, rcond=None)[0])
		return np.sqrt(
		np.clip(self.alloc_inverse @ vec, a_min=0., a_max=None)

multirotor/trajectories.py

+17 −4

Changes for multirotor/trajectories.py: 17 added lines, 4 removed lines.
Original line number	Diff line number	Diff line
		@@ -48,7 +48,11 @@ class Trajectory:
		create intermediate points a distance 2 apart.
		"""
		self.vehicle = vehicle
		self.points = np.asarray(points, self.vehicle.dtype)
		if self.vehicle is not None:
		dtype = self.vehicle.dtype
		else:
		dtype = np.float32
		self.points = np.asarray(points, dtype)
		self.proximity = proximity
		self.resolution = resolution
		self.ref = None
		@@ -95,9 +99,18 @@ class Trajectory:
		if self.resolution is not None:
		_points = []
		for i, (p1, p2) in enumerate(zip(points[:-1], points[1:])):
		dist = np.linalg.norm(p2 - p1)
		num = int(dist / self.resolution) + 1
		_points.extend(np.linspace(p1, p2, num=num, endpoint=True))
		pos_vec = p2 - p1
		dist = np.linalg.norm(pos_vec)
		unit_vec = pos_vec / dist
		# num = int(dist / self.resolution) + 1
		number = dist // self.resolution
		remainder = dist % self.resolution
		pts = p1 + unit_vec * self.resolution * np.arange(0, number+1, step=1).reshape(-1,1)
		if remainder > 0 and i==len(points)-2: # if leftover distance and last point
		pts = np.concatenate((pts, (p2,)))
		# pts = p1 + unit_vec * np.arange(0, dist + self.resolution, step=self.resolution).reshape(-1,1)
		# _points.extend(np.linspace(p1, p2, num=num, endpoint=True))
		_points.extend(pts)
		else:
		_points = points
		return _points

multirotor/vehicle.py

+2 −34

File changed.

Preview size limit exceeded, changes collapsed.

Admin message