# RCAIDE/Library/Missions/Segments/converge.py
#
#
# Created: Jul 2023, M. Clarke
# ----------------------------------------------------------------------------------------------------------------------
# IMPORT
# ----------------------------------------------------------------------------------------------------------------------
import RCAIDE
from RCAIDE.Framework.Core import Units, Data
from RCAIDE.Framework.Optimization.Packages.scipy import scipy_setup
from RCAIDE.Framework.Optimization.Common import Nexus
from RCAIDE.Framework.Analyses.Process import Process
import scipy
import scipy.optimize
import numpy as np
import sys
import os
# ----------------------------------------------------------------------------------------------------------------------
# converge root
# ----------------------------------------------------------------------------------------------------------------------
[docs]
def converge(segment):
"""Interfaces the mission a root finder algorithm.
Assumptions:
N/A
Source:
N/A
Inputs:
segment [Data]
segment.settings.root_finder [Data]
state.numerics.tolerance_solution [Unitless]
Outputs:
state.unknowns [Any]
segment.state.numerics.converged [Unitless]
Properties Used:
N/A
"""
if segment.state.numerics.solver.type == "optimize":
problem = add_mission_variables(segment)
# Commense suppression of console window output
devnull = open(os.devnull,'w')
sys.stdout = devnull
outputs = scipy_setup.SciPy_Solve(problem,
solver = segment.state.numerics.solver.method,
sense_step = segment.state.numerics.solver.step_size,
iter = segment.state.numerics.solver.max_evaluations,
tolerance = segment.state.numerics.solver.tolerance_solution)
# Terminate suppression of console window output
sys.stdout = sys.__stdout__
if outputs[3] != 0:
mission_converge = False
error_message = outputs[4]
else:
mission_converge = True
elif segment.state.numerics.solver.type == "root_finder":
unknowns = segment.state.unknowns.pack_array()
if segment.state.number_of_unknowns != segment.state.number_of_residuals:
raise AttributeError('\n The system of equations representing the mission is not square. The number of unknowns (' + str(segment.state.number_of_unknowns) + \
') is not equal to the number of residuals (equations) (' + str(segment.state.number_of_residuals) + '). Either enforce of unknowns '+\
' to be equal to the number of residuals (equations) to use fsolve or switch RCAIDE solver type to "optimize" when defining the segment.'+ \
'\n i.e. segment.state.numerics.solver.type = "optimize" ')
else:
unknowns,infodict,ier,error_message = scipy.optimize.fsolve(iterate_root_finder,
unknowns,
args = segment,
xtol = segment.state.numerics.solver.tolerance_solution,
maxfev = segment.state.numerics.solver.max_evaluations,
epsfcn = segment.state.numerics.solver.step_size,
full_output = 1)
if ier !=1:
mission_converge = False
else:
mission_converge = True
else:
raise Exception('undefined mission solver type')
if mission_converge == False:
print("Segment did not converge. Segment Tag: " + segment.tag)
print("Error Message:\n" + error_message)
segment.state.numerics.solver.converged = False
segment.converged = False
else:
segment.state.numerics.solver.converged = True
segment.converged = True
return
# ----------------------------------------------------------------------------------------------------------------------
# Helper Functions
# ----------------------------------------------------------------------------------------------------------------------
[docs]
def iterate_root_finder(unknowns, segment):
"""Runs one iteration of of all analyses for the mission.
Assumptions:
N/A
Source:
N/A
Inputs:
state.unknowns [Data]
segment.process.iterate [Data]
Outputs:
residuals [Unitless]
Properties Used:
N/A
"""
if isinstance(unknowns,np.ndarray):
segment.state.unknowns.unpack_array(unknowns)
else:
segment.state.unknowns = unknowns
segment.process.iterate(segment)
residuals = segment.state.residuals.pack_array()
return residuals
[docs]
def add_mission_variables(segment):
"""Make a pretty table view of the problem with objective and constraints at the current inputs for the dummy solver
Assumptions:
N/A
Source:
N/A
Inputs:
x [vector]
Outputs:
input [array]
const_table [array]
Properties Used:
None
"""
# Step 1: Define Nexus
nexus = Nexus()
optimization_problem = Data()
# Step 2 : Get segment type
ground_seg_flag = (type(segment) == RCAIDE.Framework.Mission.Segments.Ground.Landing) or\
(type(segment) == RCAIDE.Framework.Mission.Segments.Ground.Takeoff) or \
(type(segment) == RCAIDE.Framework.Mission.Segments.Ground.Ground)
single_pt_seg = (type(segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Altitude) or\
(type(segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Altitude_AVL_Trimmed) or \
(type(segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Altitude_No_Propulsion) or \
(type(segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Throttle)
# Step 2: Optimizer Inputs
# Step 2.1: Extract inputs
input_count = 0
unknown_keys = list(segment.state.unknowns.keys())
unknown_keys.remove('tag')
if ground_seg_flag:
n_points = segment.state.numerics.number_of_control_points
len_inputs = n_points
len_residuals = n_points
elif single_pt_seg:
n_points = 1
len_inputs = segment.state.number_of_unknowns
len_residuals = segment.state.number_of_residuals
else:
n_points = segment.state.numerics.number_of_control_points
len_inputs = n_points*segment.state.number_of_unknowns
len_residuals = n_points*segment.state.number_of_residuals
full_unkn_vals = Data()
full_upper_bound_vals = Data()
full_lower_bound_vals = Data()
for unkn in unknown_keys:
full_unkn_vals[unkn] = segment.state.unknowns[unkn]
full_lower_bound_vals[unkn] = np.atleast_2d(segment.state.numerics.solver.lower_bounds[unkn])
full_upper_bound_vals[unkn] = np.atleast_2d(segment.state.numerics.solver.upper_bounds[unkn])
# Step 2.2: Construct nexus format : [Variable_###, initial, -np.inf, np.inf , scaling, Units.less]
initial_values = full_unkn_vals.pack_array()
input_len_strings = np.tile('Variable_', len_inputs)
input_numbers = np.linspace(1,len_inputs,len_inputs,dtype=np.int16)
input_names = np.char.add(input_len_strings,np.array(input_numbers+input_count).astype(str))
lower_bounds = full_lower_bound_vals.pack_array()
upper_bounds = full_upper_bound_vals.pack_array()
units = np.broadcast_to(Units.less,(len_inputs,))
new_inputs = np.reshape(np.tile(np.atleast_2d(np.array([None,None,None,None,None,None])),len_inputs), (-1, 6))
# scaling factor for optimizer
factor = np.ceil(np.log10(abs(initial_values)))
factor[np.isinf(factor)] = 0
scale = 10 ** (factor)
# Step 2.4 Add in the inputs
new_inputs[:,0] = input_names
new_inputs[:,1] = initial_values
new_inputs[:,2] = lower_bounds
new_inputs[:,3] = upper_bounds
new_inputs[:,4] = scale
new_inputs[:,5] = units
optimization_problem.inputs = np.array(new_inputs,dtype=object)
# Step 3: Constraints
# Step 3.1 : Create the equality constraints to the beginning of the constraints all equality constraints are 0, scale 1, and unitless
new_con = np.reshape(np.tile(np.atleast_2d(np.array([None,None,None,None,None])),len_residuals), (-1, 5))
con_count = 0
con_len_strings = np.tile('Residual_', len_residuals)
con_numbers = np.linspace(1,len_residuals,len_residuals,dtype=np.int16)
con_names = np.char.add(con_len_strings,np.array(con_numbers+con_count).astype(str))
equals = np.broadcast_to('=',(len_residuals,))
zeros = np.zeros(len_residuals)
ones = np.ones(len_residuals)
# Step 3.2 Add in the new constraints
new_con[:,0] = con_names
new_con[:,1] = equals
new_con[:,2] = zeros
new_con[:,3] = ones
new_con[:,4] = 1*Units.less
optimization_problem.constraints = np.array(new_con,dtype=object)
# Step 4. Aliases
# Step 4.1: Setup the aliases for the inputs
basic_string_con = Data()
input_string = []
if ground_seg_flag:
output_numbers = np.linspace(0,n_points-2,n_points-1,dtype=np.int16)
basic_string_con[unknown_keys[1]] = np.tile('segment.state.unknowns.'+unknown_keys[1]+'[', n_points-1)
input_string.append(np.char.add(basic_string_con[unknown_keys[1]],np.array(output_numbers).astype(str)))
input_string = np.array(input_string[0])
input_string = np.char.add(input_string, np.tile(']',len_inputs-1))
input_aliases = np.reshape(np.tile(np.atleast_2d(np.array((None,None))),len_inputs), (-1, 2))
input_aliases[:,0] = input_names
input_aliases[0,1] = 'segment.state.unknowns.'+unknown_keys[0]
input_aliases[1:,1] = input_string
elif single_pt_seg:
input_string = np.array([f'segment.state.unknowns.{unkn}[0]'
for unkn in unknown_keys], dtype=str)
input_aliases = np.reshape(np.tile(np.atleast_2d(np.array((None,None))),len_inputs), (-1, 2))
input_aliases[:,0] = input_names
input_aliases[:,1] = input_string
else:
input_string = np.array([f'segment.state.unknowns.{unkn}[{i}]'
for unkn in unknown_keys
for i in range(n_points)], dtype=str)
input_aliases = np.reshape(np.tile(np.atleast_2d(np.array((None,None))),len_inputs), (-1, 2))
input_aliases[:,0] = input_names
input_aliases[:,1] = input_string
# Step 4.2: Setup the aliases for the residuals
basic_string_res = np.tile('segment.state.residuals.pack_array()[', len_residuals)
residual_string = np.char.add(basic_string_res,np.array(con_numbers-1).astype(str))
residual_string = np.char.add(residual_string, np.tile(']',len_residuals))
residual_aliases = np.reshape(np.tile(np.atleast_2d(np.array((None,None))),len_residuals), (-1, 2))
residual_aliases[:,0] = con_names
residual_aliases[:,1] = residual_string
# Step 4.3: Append Aliases
aliases = []
for ii in range(len_inputs):
aliases.append(input_aliases[ii].tolist())
for jj in range(len_residuals):
aliases.append(residual_aliases[jj].tolist())
# Step 5: Objective function
if segment.state.numerics.solver.objective == None:
aliases.append([ 'nothing' , 'postprocess.nothing'])
optimization_problem.objective = np.array([ [ 'nothing' , 1 , 1*Units.less] ],dtype=object)
elif segment.state.numerics.solver.objective == "energy":
aliases.append([ 'energy_consumed' , 'postprocess.energy_consumed'])
optimization_problem.objective = np.array([ [ 'energy_consumed' , 1 , 1*Units.less] ],dtype=object)
elif segment.state.numerics.solver.objective == "power":
aliases.append([ 'maximum_power' , 'postprocess.maximum_power'])
optimization_problem.objective = np.array([ [ 'maximum_power' , 1 , 1*Units.less] ],dtype=object)
else:
raise Exception('undefined objective function')
# append aliases
optimization_problem.aliases = aliases
# Step 6: Expand Rows
segment.process.initialize.expand_state(segment)
# Step 7: Update iteration
input_count = input_count+input_numbers[-1]
# Step 8: Append segment
nexus.segment = segment
# Step 9: Append procedure
nexus.procedure = iterate_segment()
# Step 9: Append post-process
nexus.postprocess = Data()
# Step 10: Append optimization problem
nexus.optimization_problem = optimization_problem
return nexus
[docs]
def iterate_segment():
procedure = Process()
procedure.segment = Process()
procedure.segment.design_mission = iterate_optimizer
procedure.post_process = segment_post_process
return procedure
[docs]
def iterate_optimizer(nexus):
segment = nexus.segment
unknowns = segment.state.unknowns.pack_array()
if isinstance(unknowns,np.ndarray):
segment.state.unknowns.unpack_array(unknowns)
else:
segment.state.unknowns = unknowns
segment.process.iterate(segment)
residuals = segment.state.residuals.pack_array()
nexus.residuals = residuals
return nexus
[docs]
def segment_post_process(nexus):
# unpack
power = nexus.segment.state.conditions.energy.power
I = nexus.segment.state.numerics.time.integrate
# compute max power of segment
max_power = np.max(nexus.segment.state.conditions.energy.power)
# compute total energy consumed
if (type(nexus.segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Altitude) or\
(type(nexus.segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Altitude_AVL_Trimmed) or \
(type(nexus.segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Altitude_No_Propulsion) or \
(type(nexus.segment) == RCAIDE.Framework.Mission.Segments.Single_Point.Set_Speed_Set_Throttle):
energy_consumed = 0
else:
energy_consumed = np.dot(I,power)[-1][0]
postprocess = nexus.postprocess
postprocess.maximum_power = max_power
postprocess.energy_consumed = energy_consumed
postprocess.nothing = 0
return nexus
