# ISC License
#
# Copyright (c) 2016, 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.
#
# Unit Test Script
# Module Name: cssWlsEst()
# Author: Hanspeter Schaub
# Creation Date: April 29, 2018
#
import inspect
import logging
import math
import os
import matplotlib.pyplot as plt
import numpy
import pytest
from Basilisk.architecture import messaging
from Basilisk.fswAlgorithms import cssWlsEst
# Import all of the modules that we are going to be called in this simulation
from Basilisk.utilities import SimulationBaseClass
from Basilisk.utilities import macros
from Basilisk.utilities import unitTestSupport # general support file with common unit test functions
filename = inspect.getframeinfo(inspect.currentframe()).filename
path = os.path.dirname(os.path.abspath(filename))
# Function that takes a sun pointing vector and array of CSS normal vectors and
# returns the measurements associated with those normal vectors.
def createCosList(sunPointVec, sensorPointList):
outList = []
for sensorPoint in sensorPointList:
outList.append(numpy.dot(sunPointVec, sensorPoint))
if(outList[-1] < 0.0):
outList[-1] = 0.0
return outList
# Method that checks that all of the numActive outputs from the data array
# that are greater than threshold thresh are consistent with the values in
# measVec
def checkNumActiveAccuracy(measVec, numActiveUse, numActiveFailCriteria, thresh):
numActivePred = 0
testFailCount = 0
# Iterate through measVec and find all valid signals
for i in range(0, 32):
obsVal = measVec.CosValue[i]
if (obsVal > thresh):
numActivePred += 1
# Iterate through the numActive array and sum up all numActive estimates
numActiveTotal = numpy.array([0.])
j = 0
while j < numActiveUse.shape[0]:
numActiveTotal += numActiveUse[j, 1:]
j += 1
numActiveTotal /= j # Mean number of numActive
# If we violate the test criteria, increment the failure count and alert user
if (abs(numActiveTotal[0] - numActivePred) > numActiveFailCriteria):
testFailCount += 1
errorString = "Active number failure for count of: "
errorString += str(numActivePred)
logging.error(errorString)
return testFailCount
# This method takes the sHat estimate output by the estimator and compares that
# against the actual sun vector passed in as an argument. If it doesn't match
# to the specified tolerance, increment failure counter and alert the user
def checksHatAccuracy(testVec, sHatEstUse, angleFailCriteria, TotalSim):
j = 0
testFailCount = 0
sHatTotal = numpy.array([0.0, 0.0, 0.0])
# Sum up all of the sHat estimates from the execution
while j < sHatEstUse.shape[0]:
sHatTotal += sHatEstUse[j]
j += 1
sHatTotal /= j # mean sHat estimate
# This logic is to protect cases where the dot product numerically breaks acos
dot_value = numpy.dot(sHatTotal, testVec)
if (abs(dot_value > 1.0)):
dot_value -= 2.0 * (dot_value - math.copysign(1.0, dot_value))
# If we violate the failure criteria, increment failure count and alert user
if (abs(math.acos(dot_value)) > angleFailCriteria):
testFailCount += 1
errorString = "Angle fail criteria violated for test vector:"
errorString += str(testVec).strip('[]') + "\n"
errorString += "Criteria violation of: "
errorString += str(abs(math.acos(numpy.dot(sHatTotal, testVec))))
logging.error(errorString)
return testFailCount
# This method takes the sHat estimate output by the estimator and compares that
# against the actual sun vector passed in as an argument. If it doesn't match
# to the specified tolerance, increment failure counter and alert the user
def checkResidAccuracy(testVec, sResids, sThresh, TotalSim):
j = 0
testFailCount = 0
# Sum up all of the sHat estimates from the execution
while j < sResids.shape[0]:
sNormObs = numpy.linalg.norm(sResids[j,1:])
if(sNormObs > sThresh):
testFailCount += 1
errorString = "Residual error computation failure:"
errorString += str(testVec).strip('[]') + "\n"
errorString += "Criteria violation of: "
errorString += str(sNormObs)
logging.error(errorString)
j += 1
return testFailCount
# uncomment this line is this test is to be skipped in the global unit test run, adjust message as needed
# @pytest.mark.skipif(conditionstring)
# uncomment this line if this test has an expected failure, adjust message as needed
# @pytest.mark.xfail(conditionstring)
[docs]
@pytest.mark.parametrize("testSunHeading, testRate", [
("True", "False")
,("False", "True")
])
# provide a unique test method name, starting with test_
def test_module(show_plots, testSunHeading, testRate): # update "module" in this function name to reflect the module name
"""Module Unit Test"""
# each test method requires a single assert method to be called
# pass on the testPlotFixture so that the main test function may set the DataStore attributes
if testSunHeading:
[testResults, testMessage] = cssWlsEstTestFunction(show_plots)
assert testResults < 1, testMessage
if testRate:
[testResults, testMessage] = cssRateTestFunction(show_plots)
assert testResults < 1, testMessage
def cssWlsEstTestFunction(show_plots):
testFailCount = 0 # zero unit test result counter
testMessages = [] # create empty array to store test log messages
unitTaskName = "unitTask" # arbitrary name (don't change)
unitProcessName = "TestProcess" # arbitrary name (don't change)
# Create a sim module as an empty container
unitTestSim = SimulationBaseClass.SimBaseClass()
# Create test thread
testProc = unitTestSim.CreateNewProcess(unitProcessName)
testProc.addTask(unitTestSim.CreateNewTask(unitTaskName, int(1E8)))
# Construct algorithm and associated C++ container
CSSWlsEstFSW = cssWlsEst.cssWlsEst()
CSSWlsEstFSW.ModelTag = "CSSWlsEst"
# Add module to runtime call list
unitTestSim.AddModelToTask(unitTaskName, CSSWlsEstFSW)
# Initialize the WLS estimator configuration data
CSSWlsEstFSW.useWeights = False
CSSWlsEstFSW.sensorUseThresh = 0.15
CSSOrientationList = [
[0.70710678118654746, -0.5, 0.5],
[0.70710678118654746, -0.5, -0.5],
[0.70710678118654746, 0.5, -0.5],
[0.70710678118654746, 0.5, 0.5],
[-0.70710678118654746, 0, 0.70710678118654757],
[-0.70710678118654746, 0.70710678118654757, 0.0],
[-0.70710678118654746, 0, -0.70710678118654757],
[-0.70710678118654746, -0.70710678118654757, 0.0],
]
numCSS = len(CSSOrientationList)
# set the CSS unit vectors
cssConfigData = messaging.CSSConfigMsgPayload()
totalCSSList = []
for CSSHat in CSSOrientationList:
CSSConfigElement = messaging.CSSUnitConfigMsgPayload()
CSSConfigElement.CBias = 1.0
CSSConfigElement.nHat_B = CSSHat
totalCSSList.append(CSSConfigElement)
cssConfigData.nCSS = numCSS
cssConfigData.cssVals = totalCSSList
cssConfigDataInMsg = messaging.CSSConfigMsg().write(cssConfigData)
# Initialize CSS input message
cssDataMsg = messaging.CSSArraySensorMsgPayload()
cssDataInMsg = messaging.CSSArraySensorMsg().write(cssDataMsg)
angleFailCriteria = 17.5 * math.pi / 180.0 # Get 95% effective charging in this case
numActiveFailCriteria = 0.000001 # basically zero
residFailCriteria = 1.0E-12 # Essentially numerically "small"
# Log the output message as well as the internal numACtiveCss variables
navData = CSSWlsEstFSW.navStateOutMsg.recorder()
filterData = CSSWlsEstFSW.cssWLSFiltResOutMsg.recorder()
numActiveData = CSSWlsEstFSW.logger("numActiveCss")
unitTestSim.AddModelToTask(unitTaskName, navData)
unitTestSim.AddModelToTask(unitTaskName, filterData)
unitTestSim.AddModelToTask(unitTaskName, numActiveData)
# connect the messages
CSSWlsEstFSW.cssDataInMsg.subscribeTo(cssDataInMsg)
CSSWlsEstFSW.cssConfigInMsg.subscribeTo(cssConfigDataInMsg)
# Initial test is all of the principal body axes
TestVectors = [[-1.0, 0.0, 0.0],
[0.0, -1.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, -1.0],
[0.0, 0.0, 1.0]]
# Initialize test and then step through all of the test vectors in a loop
unitTestSim.InitializeSimulation()
CSSWlsEstFSW.Reset(0) # this module reset function needs a time input (in NanoSeconds)
stepCount = 0
logLengthPrev = 0
truthData = []
for testVec in TestVectors:
if stepCount > 1: # Doing this to test permutations and get code coverage
CSSWlsEstFSW.useWeights = True
# Get observation data based on sun pointing and CSS orientation data
cssDataMsg.CosValue = createCosList(testVec, CSSOrientationList)
# Write in the observation data to the input message
cssDataInMsg.write(cssDataMsg)
# Increment the stop time to new termination value
unitTestSim.ConfigureStopTime(int((stepCount + 1) * 1E9))
# Execute simulation to current stop time
unitTestSim.ExecuteSimulation()
stepCount += 1
# Pull logged data out into workspace for analysis
sHatEst = navData.vehSunPntBdy
numActive = unitTestSupport.addTimeColumn(numActiveData.times(), numActiveData.numActiveCss)
sHatEstUse = sHatEst[logLengthPrev:, :] # Only data for this subtest
numActiveUse = numActive[logLengthPrev + 1:, :] # Only data for this subtest
# Check failure criteria and add test failures
testFailCount += checksHatAccuracy(testVec, sHatEstUse, angleFailCriteria,
unitTestSim)
testFailCount += checkNumActiveAccuracy(cssDataMsg, numActiveUse,
numActiveFailCriteria, CSSWlsEstFSW.sensorUseThresh)
filtRes = filterData.postFitRes
testFailCount += checkResidAccuracy(testVec, filtRes, residFailCriteria,
unitTestSim)
# Pop truth state onto end of array for plotting purposes
currentRow = [sHatEstUse[0, 0]]
currentRow.extend(testVec)
truthData.append(currentRow)
currentRow = [sHatEstUse[-1, 0]]
currentRow.extend(testVec)
truthData.append(currentRow)
logLengthPrev = sHatEst.shape[0]
# Hand construct case where we get low coverage (2 valid sensors)
LonVal = 0.0
LatVal = 40.68 * math.pi / 180.0
doubleTestVec = [math.sin(LatVal), math.cos(LatVal) * math.sin(LonVal),
math.cos(LatVal) * math.cos(LonVal)]
cssDataMsg.CosValue = createCosList(doubleTestVec, CSSOrientationList)
# Write in double coverage conditions and ensure that we get correct outputs
cssDataInMsg.write(cssDataMsg)
unitTestSim.ConfigureStopTime(int((stepCount + 1) * 1E9))
unitTestSim.ExecuteSimulation()
stepCount += 1
sHatEst = navData.vehSunPntBdy
numActive = unitTestSupport.addTimeColumn(numActiveData.times(), numActiveData.numActiveCss)
sHatEstUse = sHatEst[logLengthPrev:, :]
numActiveUse = numActive[logLengthPrev + 1:, :]
logLengthPrev = sHatEst.shape[0]
currentRow = [sHatEstUse[0, 0]]
currentRow.extend(doubleTestVec)
truthData.append(currentRow)
currentRow = [sHatEstUse[-1, 0]]
currentRow.extend(doubleTestVec)
truthData.append(currentRow)
# Check test criteria again
testFailCount += checksHatAccuracy(doubleTestVec, sHatEstUse, angleFailCriteria,
unitTestSim)
testFailCount += checkNumActiveAccuracy(cssDataMsg, numActiveUse,
numActiveFailCriteria, CSSWlsEstFSW.sensorUseThresh)
# Same test as above, but zero first element to get to a single coverage case
cssDataMsg.CosValue[0] = 0.0
cssDataInMsg.write(cssDataMsg)
unitTestSim.ConfigureStopTime(int((stepCount + 1) * 1E9))
unitTestSim.ExecuteSimulation()
stepCount += 1
numActive = unitTestSupport.addTimeColumn(numActiveData.times(), numActiveData.numActiveCss)
numActiveUse = numActive[logLengthPrev + 1:, :]
sHatEst = navData.vehSunPntBdy
sHatEstUse = sHatEst[logLengthPrev + 1:, :]
logLengthPrev = sHatEst.shape[0]
testFailCount += checkNumActiveAccuracy(cssDataMsg, numActiveUse,
numActiveFailCriteria, CSSWlsEstFSW.sensorUseThresh)
currentRow = [sHatEstUse[0, 0]]
currentRow.extend(doubleTestVec)
truthData.append(currentRow)
currentRow = [sHatEstUse[-1, 0]]
currentRow.extend(doubleTestVec)
truthData.append(currentRow)
# Same test as above, but zero first and fourth elements to get to zero coverage
cssDataMsg.CosValue[0] = 0.0
cssDataMsg.CosValue[3] = 0.0
cssDataInMsg.write(cssDataMsg)
unitTestSim.ConfigureStopTime(int((stepCount + 1) * 1E9))
unitTestSim.ExecuteSimulation()
numActive = unitTestSupport.addTimeColumn(numActiveData.times(), numActiveData.numActiveCss)
numActiveUse = numActive[logLengthPrev:, :]
logLengthPrev = numActive.shape[0]
testFailCount += checkNumActiveAccuracy(cssDataMsg, numActiveUse,
numActiveFailCriteria, CSSWlsEstFSW.sensorUseThresh)
# Format data for plotting
truthData = numpy.array(truthData)
sHatEst = navData.vehSunPntBdy
numActive = unitTestSupport.addTimeColumn(numActiveData.times(), numActiveData.numActiveCss)
#
# test the case where all CSS signals are zero
#
cssDataMsg.CosValue = numpy.zeros(len(CSSOrientationList))
cssDataInMsg.write(cssDataMsg)
unitTestSim.ConfigureStopTime(int((stepCount + 2) * 1E9))
unitTestSim.ExecuteSimulation()
sHatEstZero = navData.vehSunPntBdy
sHatEstZeroUse = sHatEstZero[logLengthPrev + 1:, :]
trueVector = [[0.0, 0.0, 0.0]]*len(sHatEstZeroUse)
for i in range(0, len(trueVector)):
# check a vector values
if not unitTestSupport.isArrayEqual(sHatEstZeroUse[i], trueVector[i], 3, 1e-12):
testFailCount += 1
testMessages.append("FAILED: " + CSSWlsEstFSW.ModelTag + " Module failed unit test at t=" +
str(navData.times()[i] * macros.NANO2SEC) + "sec\n")
if show_plots:
plt.figure(1)
plt.plot(sHatEst[:, 0] * 1.0E-9, sHatEst[:, 0], label='x-Sun')
plt.plot(sHatEst[:, 0] * 1.0E-9, sHatEst[:, 1], label='y-Sun')
plt.plot(sHatEst[:, 0] * 1.0E-9, sHatEst[:, 2], label='z-Sun')
plt.legend(loc='upper left')
plt.xlabel('Time (s)')
plt.ylabel('Unit Component (--)')
plt.figure(2)
plt.plot(numActive[:, 0] * 1.0E-9, numActive[:, 1])
plt.xlabel('Time (s)')
plt.ylabel('Number Active CSS (--)')
plt.figure(3)
plt.subplot(3, 1, 1)
plt.plot(sHatEst[:, 0] * 1.0E-9, sHatEst[:, 0], label='Est')
plt.plot(truthData[:, 0] * 1.0E-9, truthData[:, 0], 'r--', label='Truth')
plt.xlabel('Time (s)')
plt.ylabel('X Component (--)')
plt.legend(loc='lower right')
plt.subplot(3, 1, 2)
plt.plot(sHatEst[:, 0] * 1.0E-9, sHatEst[:, 1], label='Est')
plt.plot(truthData[:, 0] * 1.0E-9, truthData[:, 1], 'r--', label='Truth')
plt.xlabel('Time (s)')
plt.ylabel('Y Component (--)')
plt.subplot(3, 1, 3)
plt.plot(sHatEst[:, 0] * 1.0E-9, sHatEst[:, 2], label='Est')
plt.plot(truthData[:, 0] * 1.0E-9, truthData[:, 2], 'r--', label='Truth')
plt.xlabel('Time (s)')
plt.ylabel('Z Component (--)')
plt.show()
plt.close('all')
# print out success message if no error were found
if testFailCount == 0:
print("PASSED: " + CSSWlsEstFSW.ModelTag)
# each test method requires a single assert method to be called
# this check below just makes sure no sub-test failures were found
return [testFailCount, ''.join(testMessages)]
def cssRateTestFunction(show_plots):
testFailCount = 0 # zero unit test result counter
testMessages = [] # create empty array to store test log messages
unitTaskName = "unitTask" # arbitrary name (don't change)
unitProcessName = "TestProcess" # arbitrary name (don't change)
# Create a sim module as an empty container
unitTestSim = SimulationBaseClass.SimBaseClass()
# Create test thread
testProc = unitTestSim.CreateNewProcess(unitProcessName)
testProcessRate = macros.sec2nano(0.5) # update process rate update time
testProc.addTask(unitTestSim.CreateNewTask(unitTaskName, testProcessRate))
# Construct algorithm and associated C++ container
module = cssWlsEst.cssWlsEst()
module.ModelTag = "CSSWlsEst"
# Add module to runtime call list
unitTestSim.AddModelToTask(unitTaskName, module)
# Initialize the WLS estimator configuration data
module.useWeights = False
module.sensorUseThresh = 0.15
CSSOrientationList = [
[0.70710678118654746, -0.5, 0.5],
[0.70710678118654746, -0.5, -0.5],
[0.70710678118654746, 0.5, -0.5],
[0.70710678118654746, 0.5, 0.5],
[-0.70710678118654746, 0, 0.70710678118654757],
[-0.70710678118654746, 0.70710678118654757, 0.0],
[-0.70710678118654746, 0, -0.70710678118654757],
[-0.70710678118654746, -0.70710678118654757, 0.0],
]
numCSS = len(CSSOrientationList)
# set the CSS unit vectors
cssConfigData = messaging.CSSConfigMsgPayload()
totalCSSList = []
for CSSHat in CSSOrientationList:
CSSConfigElement = messaging.CSSUnitConfigMsgPayload()
CSSConfigElement.CBias = 1.0
CSSConfigElement.nHat_B = CSSHat
totalCSSList.append(CSSConfigElement)
cssConfigData.nCSS = numCSS
cssConfigData.cssVals = totalCSSList
cssConfigDataInMsg = messaging.CSSConfigMsg().write(cssConfigData)
# Log the output message as well as the internal numACtiveCss variables
dataLog = module.navStateOutMsg.recorder()
unitTestSim.AddModelToTask(unitTaskName, dataLog)
# Get observation data based on sun pointing and CSS orientation data
cssDataMsg = messaging.CSSArraySensorMsgPayload()
cssDataMsg.CosValue = createCosList([1.0, 0.0, 0.0], CSSOrientationList)
cssDataInMsg = messaging.CSSArraySensorMsg().write(cssDataMsg)
# connect messages
module.cssDataInMsg.subscribeTo(cssDataInMsg)
module.cssConfigInMsg.subscribeTo(cssConfigDataInMsg)
# Initialize test and then step through all of the test vectors in a loop
unitTestSim.InitializeSimulation()
# Increment the stop time to new termination value
unitTestSim.ConfigureStopTime(macros.sec2nano(1.0))
# Execute simulation to current stop time
unitTestSim.ExecuteSimulation()
# rotate sun heading by 90 degrees
cssDataMsg.CosValue = createCosList([0.0, 1.0, 0.0], CSSOrientationList)
cssDataInMsg.write(cssDataMsg)
unitTestSim.ConfigureStopTime(macros.sec2nano(2.0))
unitTestSim.ExecuteSimulation()
# test the module reset function
module.Reset(1) # this module reset function needs a time input (in NanoSeconds)
unitTestSim.ConfigureStopTime(macros.sec2nano(2.5))
unitTestSim.ExecuteSimulation()
cssDataMsg.CosValue = createCosList([1.0, 0.0, 0.0], CSSOrientationList)
cssDataInMsg.write(cssDataMsg)
unitTestSim.ConfigureStopTime(macros.sec2nano(3.0))
unitTestSim.ExecuteSimulation()
accuracy = 1e-6
trueVector = [
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, -3.14159265],
[0.0, 0.0, 0.0],
[0.0, 0.0, 0.0],
[0.0, 0.0, +3.14159265]
]
testFailCount, testMessages = unitTestSupport.compareArray(trueVector, dataLog.omega_BN_B,
accuracy, "CSS Rate Vector",
testFailCount, testMessages)
# print out success message if no error were found
snippentName = "passFailRate"
if testFailCount == 0:
colorText = 'ForestGreen'
print("PASSED: " + module.ModelTag)
passedText = r'\textcolor{' + colorText + '}{' + "PASSED" + '}'
else:
colorText = 'Red'
print("Failed: " + module.ModelTag)
passedText = r'\textcolor{' + colorText + '}{' + "Failed" + '}'
unitTestSupport.writeTeXSnippet(snippentName, passedText, path)
# each test method requires a single assert method to be called
# this check below just makes sure no sub-test failures were found
return [testFailCount, ''.join(testMessages)]
#
# This statement below ensures that the unitTestScript can be run as a
# stand-along python script
#
if __name__ == "__main__":
test_module(
True, # show_plots
False, # testSunHeading Flag
True # testRate Flag
)