#
# ISC License
#
# Copyright (c) 2022, Autonomous Vehicle Systems Lab, University of Colorado at Boulder
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
#
r"""
Overview
--------
This scenario demonstrates how to set up a spacecraft with deploying panels. The panel modules
are set up with individual profilers and motors as shown in the following diagram.
.. image:: /_images/static/test_scenario_DeployingPanel.svg
:align: center
The script is found in the folder ``basilisk/examples`` and executed by using::
python3 scenarioDeployingPanel.py
The simulation includes two deploying panels that start undeployed. The first panel deploys fully,
but the second panel deploys off-nominally (to 80%), leading to a reduced power output.
Illustration of Simulation Results
----------------------------------
::
show_plots = True
Five plots are shown. The first two show the body attitude and rate relative to the inertial frame,
demonstrating the dynamic effect of the deploying panels on the system. Because the second panel
fails to fully deploy, the body attitude does not return to the original state.
.. image:: /_images/Scenarios/scenarioDeployingPanel1.svg
:align: center
.. image:: /_images/Scenarios/scenarioDeployingPanel2.svg
:align: center
The next two plots show the panel angles and angle rates. Because the panels are modeled as
rigid bodies attached by a torsional spring, oscillations about the time-varying reference
angle are seen.
.. image:: /_images/Scenarios/scenarioDeployingPanel3.svg
:align: center
.. image:: /_images/Scenarios/scenarioDeployingPanel4.svg
:align: center
The final plot shows the power output. As the first panel deploys, power starts to be generated.
The second panel also generates power, but stops short of full deployment and thus generates less power.
.. image:: /_images/Scenarios/scenarioDeployingPanel5.svg
:align: center
"""
#
# Basilisk Scenario Script and Integrated Test
#
# Purpose: Scenario script and test of the hingedRigidBodyStateEffector with changes
# made to permit deployment.
# Author: Galen Bascom
# Creation Date: Sept. 12, 2022
#
import os
import matplotlib.pyplot as plt
import numpy as np
from Basilisk import __path__
bskPath = __path__[0]
fileName = os.path.basename(os.path.splitext(__file__)[0])
from Basilisk.simulation import spacecraft
from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
simIncludeGravBody, unitTestSupport, vizSupport)
from Basilisk.simulation import hingedRigidBodyStateEffector, simpleSolarPanel
from Basilisk.simulation import hingedBodyLinearProfiler, hingedRigidBodyMotor
import math
[docs]def run(show_plots):
"""
Args:
show_plots (bool): Determines if the script should display plots.
"""
# Create simulation variable names
simTaskName = "simTask"
simProcessName = "simProcess"
# Create a sim module as an empty container
scSim = SimulationBaseClass.SimBaseClass()
# create the simulation process
dynProcess = scSim.CreateNewProcess(simProcessName)
# create the dynamics task and specify the integration update time
simulationTimeStep = macros.sec2nano(.1)
dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))
# create the spacecraft hub
scObject = spacecraft.Spacecraft()
scObject.ModelTag = "bskSat"
scObject.hub.mHub = 750.0
scObject.hub.IHubPntBc_B = [[900.0, 0.0, 0.0], [0.0, 800.0, 0.0], [0.0, 0.0, 600.0]]
# setup Earth Gravity Body
gravFactory = simIncludeGravBody.gravBodyFactory()
gravBodies = gravFactory.createBodies('earth', 'sun')
gravBodies['earth'].isCentralBody = True
mu = gravBodies['earth'].mu
sun = 1
gravFactory.addBodiesTo(scObject)
timeInitString = "2012 MAY 1 00:28:30.0"
spiceObject = gravFactory.createSpiceInterface(time=timeInitString, epochInMsg=True)
spiceObject.zeroBase = 'earth'
scSim.AddModelToTask(simTaskName, spiceObject)
# setup the orbit using classical orbit elements
oe = orbitalMotion.ClassicElements()
rLEO = 7000. * 1000 # meters
rGEO = 42000. * 1000 # meters
oe.a = (rLEO + rGEO) / 2.0
oe.e = 1.0 - rLEO / oe.a
oe.i = 0.0 * macros.D2R
oe.Omega = 48.2 * macros.D2R
oe.omega = 347.8 * macros.D2R
oe.f = 85.3 * macros.D2R
rN, vN = orbitalMotion.elem2rv(mu, oe)
oe = orbitalMotion.rv2elem(mu, rN, vN) # this stores consistent initial orbit elements
# To set the spacecraft initial conditions, the following initial position and velocity variables are set:
scObject.hub.r_CN_NInit = rN # m - r_BN_N
scObject.hub.v_CN_NInit = vN # m/s - v_BN_N
# point the body 3 axis towards the sun in the inertial n2 direction
scObject.hub.sigma_BNInit = [[math.tan(-90. / 4. * macros.D2R)], [0.0], [0.0]] # sigma_BN_B
scObject.hub.omega_BN_BInit = [[0.0], [0.0], [0.0]] # rad/s - omega_BN_B
# configure panels
panel1 = hingedRigidBodyStateEffector.HingedRigidBodyStateEffector()
panel1.ModelTag = "panel1"
panel1.mass = 100.0
panel1.IPntS_S = [[100.0, 0.0, 0.0], [0.0, 50.0, 0.0], [0.0, 0.0, 50.0]]
panel1.d = 1.5
panel1.k = 200.
panel1.c = 20.
panel1.r_HB_B = [[-.5], [0.0], [-1.0]]
panel1.dcm_HB = [[-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0]]
panel1.thetaInit = -np.pi
panel1.thetaDotInit = 0
panel2 = hingedRigidBodyStateEffector.HingedRigidBodyStateEffector()
panel2.ModelTag = "panel2"
panel2.mass = 100.0
panel2.IPntS_S = [[100.0, 0.0, 0.0], [0.0, 50.0, 0.0], [0.0, 0.0, 50.0]]
panel2.d = 1.5
panel2.k = 200.
panel2.c = 20.
panel2.r_HB_B = [[.5], [0.0], [-1.0]]
panel2.dcm_HB = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
panel2.thetaInit = -np.pi
panel2.thetaDotInit = 0
# profilers
profiler1 = hingedBodyLinearProfiler.HingedBodyLinearProfiler()
profiler1.ModelTag = "deploymentProfiler"
profiler1.startTime = macros.sec2nano(5) # [ns] start the deployment
profiler1.endTime = macros.sec2nano(85) # [ns]
profiler1.startTheta = -np.pi # [rad] starting angle in radians
profiler1.endTheta = -np.pi / 2 # [rad]
profiler2 = hingedBodyLinearProfiler.HingedBodyLinearProfiler()
profiler2.ModelTag = "deploymentProfiler2"
profiler2.startTime = macros.sec2nano(100) # [ns] start the deployment
profiler2.endTime = macros.sec2nano(164) # [ns]
profiler2.startTheta = -np.pi # [rad] starting angle in radians
profiler2.endTheta = -np.pi / 1.75 # [rad] ending angle is not all the way deployed
panel1.hingedRigidBodyRefMsg.subscribeTo(profiler1.hingedRigidBodyReferenceOutMsg)
panel2.hingedRigidBodyRefMsg.subscribeTo(profiler2.hingedRigidBodyReferenceOutMsg)
# motors
motor1 = hingedRigidBodyMotor.HingedRigidBodyMotor()
motor1.ModelTag = "hingedRigidBodyMotor"
motor1.K = 20 # proportional gain constant
motor1.P = 10 # derivative gain constant
motor2 = hingedRigidBodyMotor.HingedRigidBodyMotor()
motor2.ModelTag = "hingedRigidBodyMotor2"
motor2.K = 20 # proportional gain constant
motor2.P = 10 # derivative gain constant
motor1.hingedBodyStateSensedInMsg.subscribeTo(panel1.hingedRigidBodyOutMsg)
motor1.hingedBodyStateReferenceInMsg.subscribeTo(profiler1.hingedRigidBodyReferenceOutMsg)
panel1.motorTorqueInMsg.subscribeTo(motor1.motorTorqueOutMsg)
motor2.hingedBodyStateSensedInMsg.subscribeTo(panel2.hingedRigidBodyOutMsg)
motor2.hingedBodyStateReferenceInMsg.subscribeTo(profiler2.hingedRigidBodyReferenceOutMsg)
panel2.motorTorqueInMsg.subscribeTo(motor2.motorTorqueOutMsg)
# add panel to spacecraft hub
scObject.addStateEffector(panel1) # in order to affect dynamics
scObject.addStateEffector(panel2) # in order to affect dynamics
# power
solarPanel1 = simpleSolarPanel.SimpleSolarPanel()
solarPanel1.ModelTag = "pwr1"
solarPanel1.nHat_B = [0, 0, 1]
solarPanel1.panelArea = 2.0 # m^2
solarPanel1.panelEfficiency = 0.9 # 90% efficiency in power generation
solarPanel1.stateInMsg.subscribeTo(panel1.hingedRigidBodyConfigLogOutMsg)
solarPanel1.sunInMsg.subscribeTo(spiceObject.planetStateOutMsgs[sun])
solarPanel2 = simpleSolarPanel.SimpleSolarPanel()
solarPanel2.ModelTag = "pwr2"
solarPanel2.nHat_B = [0, 0, 1]
solarPanel2.panelArea = 2.0 # m^2
solarPanel2.panelEfficiency = 0.9 # 90% efficiency in power generation
solarPanel2.stateInMsg.subscribeTo(panel2.hingedRigidBodyConfigLogOutMsg)
solarPanel2.sunInMsg.subscribeTo(spiceObject.planetStateOutMsgs[sun])
#
# add modules to simulation task list
#
scSim.AddModelToTask(simTaskName, scObject)
scSim.AddModelToTask(simTaskName, panel1)
scSim.AddModelToTask(simTaskName, profiler1)
scSim.AddModelToTask(simTaskName, motor1)
scSim.AddModelToTask(simTaskName, panel2)
scSim.AddModelToTask(simTaskName, profiler2)
scSim.AddModelToTask(simTaskName, motor2)
scSim.AddModelToTask(simTaskName, solarPanel1)
scSim.AddModelToTask(simTaskName, solarPanel2)
# set the simulation time
simulationTime = macros.sec2nano(240)
# Setup data logging before the simulation is initialized
numDataPoints = 1000
samplingTime = unitTestSupport.samplingTime(simulationTime, simulationTimeStep, numDataPoints)
dataLog = scObject.scStateOutMsg.recorder(samplingTime)
p1Log = panel1.hingedRigidBodyOutMsg.recorder(samplingTime)
prof1Log = profiler1.hingedRigidBodyReferenceOutMsg.recorder(samplingTime)
motor1Log = motor1.motorTorqueOutMsg.recorder(samplingTime)
p2Log = panel2.hingedRigidBodyOutMsg.recorder(samplingTime)
prof2Log = profiler2.hingedRigidBodyReferenceOutMsg.recorder(samplingTime)
motor2Log = motor2.motorTorqueOutMsg.recorder(samplingTime)
pwr1Log = solarPanel1.nodePowerOutMsg.recorder(samplingTime)
pwr2Log = solarPanel2.nodePowerOutMsg.recorder(samplingTime)
scSim.AddModelToTask(simTaskName, dataLog)
scSim.AddModelToTask(simTaskName, p1Log)
scSim.AddModelToTask(simTaskName, prof1Log)
scSim.AddModelToTask(simTaskName, motor1Log)
scSim.AddModelToTask(simTaskName, p2Log)
scSim.AddModelToTask(simTaskName, prof2Log)
scSim.AddModelToTask(simTaskName, motor2Log)
scSim.AddModelToTask(simTaskName, pwr1Log)
scSim.AddModelToTask(simTaskName, pwr2Log)
if vizSupport.vizFound:
viz = vizSupport.enableUnityVisualization(scSim, simTaskName,
[scObject
, [panel1.ModelTag, panel1.hingedRigidBodyConfigLogOutMsg]
, [panel2.ModelTag, panel2.hingedRigidBodyConfigLogOutMsg]
]
# , saveFile=__file__
)
vizSupport.createCustomModel(viz
, simBodiesToModify=[panel1.ModelTag]
, modelPath="CUBE"
, scale=[3, 1, 0.1]
, color=vizSupport.toRGBA255("blue"))
vizSupport.createCustomModel(viz
, simBodiesToModify=[panel2.ModelTag]
, modelPath="CUBE"
, scale=[3, 1, 0.1]
, color=vizSupport.toRGBA255("blue"))
viz.settings.orbitLinesOn = -1
scSim.InitializeSimulation()
scSim.ConfigureStopTime(simulationTime)
scSim.ExecuteSimulation()
# retrieve the logged data
dataSigmaBN = dataLog.sigma_BN
dataOmegaBN_B = dataLog.omega_BN_B
panel1thetaLog = p1Log.theta
panel1thetaDotLog = p1Log.thetaDot
panel2thetaLog = p2Log.theta
panel2thetaDotLog = p2Log.thetaDot
pwrLog1 = pwr1Log.netPower
pwrLog2 = pwr2Log.netPower
np.set_printoptions(precision=16)
figureList = plotOrbits(dataLog.times(), dataSigmaBN, dataOmegaBN_B,
panel1thetaLog, panel1thetaDotLog,
panel2thetaLog, panel2thetaDotLog,
pwrLog1, pwrLog2)
if show_plots:
plt.show()
# close the plots being saved off to avoid over-writing old and new figures
plt.close("all")
return figureList
def plotOrbits(timeAxis, dataSigmaBN, dataOmegaBN,
panel1thetaLog, panel1thetaDotLog,
panel2thetaLog, panel2thetaDotLog,
pwrLog1, pwrLog2):
plt.close("all") # clears out plots from earlier test runs
# sigma B/N
figCounter = 1
plt.figure(figCounter)
fig = plt.gcf()
ax = fig.gca()
ax.ticklabel_format(useOffset=False, style='plain')
timeData = timeAxis * macros.NANO2SEC
for idx in range(3):
plt.plot(timeData, dataSigmaBN[:, idx],
color=unitTestSupport.getLineColor(idx, 3),
label=r'$\sigma_' + str(idx) + '$')
plt.legend(loc='lower right')
plt.xlabel('Time [s]')
plt.ylabel(r'MRP Attitude $\sigma_{B/N}$')
figureList = {}
pltName = fileName + str(figCounter)
figureList[pltName] = plt.figure(figCounter)
# omega B/N
figCounter += 1
plt.figure(figCounter)
fig = plt.gcf()
ax = fig.gca()
ax.ticklabel_format(useOffset=False, style='plain')
for idx in range(3):
plt.plot(timeData, dataOmegaBN[:, idx],
color=unitTestSupport.getLineColor(idx, 3),
label=r'$\omega_' + str(idx) + '$')
plt.legend(loc='lower right')
plt.xlabel('Time [s]')
plt.ylabel(r'Rate $\omega_{B/N}$')
pltName = fileName + str(figCounter)
figureList[pltName] = plt.figure(figCounter)
# rotating panel hinge angle: panel 1
figCounter += 1
plt.figure(figCounter)
ax1 = plt.figure(figCounter).add_subplot(111)
ax1.plot(timeData, panel1thetaLog, color='royalblue')
plt.xlabel('Time [s]')
plt.ylabel('Panel 1 Angle [rad]', color='royalblue')
ax2 = plt.figure(figCounter).add_subplot(111, sharex=ax1, frameon=False)
ax2.plot(timeData, panel1thetaDotLog, color='indianred')
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
plt.ylabel('Panel 1 Angle Rate [rad/s]', color='indianred')
pltName = fileName + str(figCounter)
figureList[pltName] = plt.figure(figCounter)
# rotating panel hinge angle: panel 2
figCounter += 1
plt.figure(figCounter)
ax1 = plt.figure(figCounter).add_subplot(111)
ax1.plot(timeData, panel2thetaLog, color='royalblue')
plt.xlabel('Time [s]')
plt.ylabel('Panel 2 Angle [rad]', color='royalblue')
ax2 = plt.figure(figCounter).add_subplot(111, sharex=ax1, frameon=False)
ax2.plot(timeData, panel2thetaDotLog, color='indianred')
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
plt.ylabel('Panel 2 Angle Rate [rad/s]', color='indianred')
pltName = fileName + str(figCounter)
figureList[pltName] = plt.figure(figCounter)
# power: panel 1
figCounter += 1
plt.figure(figCounter)
ax1 = plt.figure(figCounter).add_subplot(111)
ax1.plot(timeData, pwrLog1, color='goldenrod', label="Panel 1")
ax1.plot(timeData, pwrLog2, '--', color='goldenrod', label="Panel 2")
plt.xlabel('Time [s]')
plt.ylabel('Panel Power [W]')
plt.legend(loc='lower right')
pltName = fileName + str(figCounter)
figureList[pltName] = plt.figure(figCounter)
return figureList
if __name__ == "__main__":
run(
True # show_plots
)