diff --git a/constraint_analysis/CMakeLists.txt b/constraint_analysis/CMakeLists.txt
index 564eb7d2292d0bdabd9d971bc4106c01dfc65462..20918bf442a22fa9ca7a5ec66ceefb867c3ad8cb 100644
--- a/constraint_analysis/CMakeLists.txt
+++ b/constraint_analysis/CMakeLists.txt
@@ -7,7 +7,7 @@ set(MODULE_SOURCES
src/constraint_analysis.cpp
src/ca_parser.cpp
src/io/aircraft_xml.cpp
- src/energy_based/energy_based.cpp
+ src/energy_based/energy_based_refactored.cpp
)
add_executable(${MODULE_NAME}
diff --git a/constraint_analysis/constraint_analysis_conf.xml b/constraint_analysis/constraint_analysis_conf.xml
index 5a7ca22c039c3af6afa6a1cd8cb437268618c4d0..bba4c8199e1a0a4a1c84a535a6e2dbe8f9c9ee16 100644
--- a/constraint_analysis/constraint_analysis_conf.xml
+++ b/constraint_analysis/constraint_analysis_conf.xml
@@ -18,7 +18,7 @@
</log_file_output>
<plot_output description="Specify the way plotting shall be handled">
<enable description="Switch to enable plotting. Switch: true (On) / false (Off)">
- <value>true</value>
+ <value>false</value>
</enable>
<copy_plotting_files description="Switch if plotting files shall be copied. Switch: true (On) / false (Off)">
<value>true</value>
@@ -50,35 +50,177 @@
</gnuplot_path>
</control_settings>
<program_settings description="Settings specific to constraint_sizing.">
- <method description="The method of constraint sizing, only Energy_Based is supported for now">
+ <method description="The method of constraint analysis, only Energy_Based is supported for now">
<value>Energy_Based</value>
</method>
- <aero_method description="The method of calculating the polar, only Read_Polar and Calculate_Polar valid.">
- <value>Calculate_Polar</value>
- </aero_method>
- <Mach_TO description="Mach number at takeoff. Caution. It has to be the Mach number number of the polar of the takeoff configuration!">
- <value>0.2</value>
- </Mach_TO>
- <takeoff_climb_angle description="Takeoff climb angle">
- <value>5</value>
- </takeoff_climb_angle>
- <gust_speed description="Gust speed at takeoff and landing normally 20 m/s">
- <value>20</value>
- </gust_speed>
- <gust_load_factor description="Additional gust load factor">
- <value>1.5</value>
- </gust_load_factor>
+ <mode_selection description="The method of sizign, mode_0 changes both T/W and W/S, mode_1 changes only the T/W">
+ <value>mode_1</value>
+ </mode_selection>
+ <buffer_factor description="The buffer in percentage, if the rate of change of the design point parameters is less than this percentage, it remains constant">
+ <value>0.0</value>
+ </buffer_factor>
<safety_factor description="Additional safety factor">
<value>0.05</value>
</safety_factor>
- <oswald_factor description="Oswald factor for calculating the polar">
- <value>0.8</value>
- </oswald_factor>
- <climb_gradient_OEI description="Minimal climb rate of 2.4% defined in CS25 for second segment climb">
- <value>0.024</value>
- </climb_gradient_OEI>
- <minimum_climb_rate description="Minimum climb rate out of CS 25: 100fpm -> 0.508 m/s">
- <value>0.508</value>
- </minimum_climb_rate>
+ <constraint_selection>
+ <standard_set description="The standard set of constraints, assembled using the requirements and certain assumptions">
+ <takeoff_ground_roll>
+ <active>
+ <value>true</value>
+ </active>
+ <Mach_TO description="Mach number at takeoff. Caution! It has to be the Mach number number of the polar of the takeoff configuration!">
+ <value>0.2</value>
+ </Mach_TO>
+ </takeoff_ground_roll>
+ <takeoff_climb>
+ <active>
+ <value>true</value>
+ </active>
+ <Mach_TO description="Mach number at takeoff. Caution! It has to be the Mach number number of the polar of the takeoff configuration!">
+ <value>0.2</value>
+ </Mach_TO>
+ <takeoff_climb_angle description="Takeoff climb angle">
+ <value>5</value>
+ </takeoff_climb_angle>
+ </takeoff_climb>
+ <OEI_climb>
+ <active>
+ <value>true</value>
+ </active>
+ <climb_gradient_OEI description="Minimal climb rate of 2.4% defined in CS25 for second segment climb">
+ <value>0.024</value>
+ </climb_gradient_OEI>
+ </OEI_climb>
+ <top_of_climb>
+ <active>
+ <value>false</value>
+ </active>
+ </top_of_climb>
+ <cruise>
+ <active>
+ <value>false</value>
+ </active>
+ </cruise>
+ <SEP>
+ <active>
+ <value>true</value>
+ </active>
+ <minimum_climb_rate description="Minimum climb rate out of CS 25: 100fpm -> 0.508 m/s">
+ <value>0.508</value>
+ </minimum_climb_rate>
+ </SEP>
+ <approach_speed>
+ <active>
+ <value>true</value>
+ </active>
+ </approach_speed>
+ <landing_field_length>
+ <active>
+ <value>true</value>
+ </active>
+ <Mach>
+ <value>0.2</value>
+ </Mach>
+ </landing_field_length>
+ <gust>
+ <active>
+ <value>true</value>
+ </active>
+ <Mach>
+ <value>0.2</value>
+ </Mach>
+ <gust_speed description="Gust speed at takeoff and landing normally 20 m/s">
+ <value>20</value>
+ </gust_speed>
+ <gust_load_factor description="Additional gust load factor">
+ <value>1.5</value>
+ </gust_load_factor>
+ </gust>
+ </standard_set>
+ <extra_set description="The extra set of constraints, assembled with user input, use the templates">
+ <cruise>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>cruise_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>10000</altitude>
+ </cruise>
+ <climb>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>climb_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>10000</altitude>
+ <climb_rate>0.508</climb_rate>
+ </climb>
+ <constant_speed_turn>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>turn_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>10000</altitude>
+ <load_factor>2</load_factor>
+ </constant_speed_turn>
+ <horizontal_acceleration>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>acceleration_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>10000</altitude>
+ <acceleration>1.0</acceleration>
+ </horizontal_acceleration>
+ <takeoff_ground_roll>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>takeoff_ground_roll_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>0</altitude>
+ <takeoff_field_length>2000</takeoff_field_length>
+ </takeoff_ground_roll>
+ <braking_roll>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>landing_field_length_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>0</altitude>
+ <friction_coefficient>0.32</friction_coefficient>
+ <landing_field_length>2000</landing_field_length>
+ </braking_roll>
+ <takeoff_climb>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>takeoff_climb_angle_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>0</altitude>
+ <takeoff_climb_angle>0.32</takeoff_climb_angle>
+ </takeoff_climb>
+ <gust>
+ <active>
+ <value>false</value>
+ </active>
+ <segment_ID>gust_template</segment_ID>
+ <weight_fraction>0.98</weight_fraction>
+ <Mach>0.8</Mach>
+ <altitude>0</altitude>
+ <V_TO_L>200</V_TO_L>
+ <w_g>20</w_g>
+ <dn_G>1.5</dn_G>
+ </gust>
+ </extra_set>
+ </constraint_selection>
</program_settings>
</module_configuration_file>
diff --git a/constraint_analysis/include/constraint_analysis/ca_functions.h b/constraint_analysis/include/constraint_analysis/ca_functions.h
index 36823d3ec4a03ee3453ce2f99e45b6233f02976b..3681889f8eee86a4aaf2110a7368c1813ae6892e 100644
--- a/constraint_analysis/include/constraint_analysis/ca_functions.h
+++ b/constraint_analysis/include/constraint_analysis/ca_functions.h
@@ -158,7 +158,7 @@ namespace Mattingly {
* @param CL_max_TO Value of the lift coefficient at takeoff.
* @param takeoff_climb_angle Value of the take-off climb angle.
*/
- auto takeoff_climb_angle(const std::vector<double>& W_over_S_data, double K_1, double K_2, double C_D0, double weight_fraction, double thrust_lapse, double Mach, double altitude, double CL_max_TO, double takeoff_climb_angle) -> const std::vector<double>;
+ auto takeoff_climb_angle(const std::vector<double>& W_over_S_data, const std::vector<double>& C_D_data, double weight_fraction, double thrust_lapse, double Mach, double altitude, double CL_max_TO, double takeoff_climb_angle) -> const std::vector<double>;
/**
* @brief Calculates and outputs thrust-to-weight ratios (T/W)
diff --git a/constraint_analysis/include/constraint_analysis/ca_minfinder.h b/constraint_analysis/include/constraint_analysis/ca_minfinder.h
index 8610d87c8515682a5da26d4ca8e9c7cce7fa474d..c4b3f8a2210fe16729482b5ad27e5afeafbeff78 100644
--- a/constraint_analysis/include/constraint_analysis/ca_minfinder.h
+++ b/constraint_analysis/include/constraint_analysis/ca_minfinder.h
@@ -20,6 +20,9 @@
* This file is part of UNICADO.
*/
+#ifndef SRC_CONSTRAINT_ANALYSIS_MINFINDER_H_
+#define SRC_CONSTRAINT_ANALYSIS_MINFINDER_H_
+
#include <cmath>
#include <vector>
#include <numbers>
@@ -107,6 +110,33 @@ public:
design_thrust_to_weight = design_thrust_to_weight * (1.0 + safety_factor);
};
+ void find_converged_design_point(const double& previous_TW, const double& previous_WS, const double& safety_factor, const double& buffer_factor)
+ {
+ std::sort(wing_loading_boundaries.begin(), wing_loading_boundaries.end());
+ std::vector<double> xlims = {};
+ for (auto x_limit : wing_loading_boundaries) {
+ auto it = std::distance(wing_loading_vec.begin(), std::lower_bound(wing_loading_vec.begin(), wing_loading_vec.end(), x_limit));
+ xlims.push_back(it);
+ }
+ std::vector<double> limited_TW(dominant_thrust_to_weight.begin() + xlims[0], dominant_thrust_to_weight.begin() + xlims[1]);
+ int i = xlims[0];
+ for (double point : limited_TW)
+ {
+ if (point <= design_thrust_to_weight)
+ {
+ design_thrust_to_weight = point;
+ design_wing_loading = wing_loading_vec[i];
+ }
+ i++;
+ }
+ design_thrust_to_weight = design_thrust_to_weight * (1.0 + safety_factor);
+
+ double relative_change_TW = sqrt(pow(1 - previous_TW / design_thrust_to_weight, 2));
+ if (relative_change_TW <= buffer_factor/100.0) { design_thrust_to_weight = previous_TW; };
+ double relative_change_WS = sqrt(pow(1 - previous_WS / design_wing_loading, 2));
+ if (relative_change_WS <= buffer_factor/100.0) { design_wing_loading = previous_WS; };
+ };
+
template<typename T>
auto create_endnode(const T value, const std::string& unit, const std::string& node_path, const std::string& description) -> Endnode<T> {
/* Create the node an set its properties */
@@ -162,3 +192,5 @@ public:
infile.close();
};
};
+
+#endif
diff --git a/constraint_analysis/include/constraint_analysis/ca_parser.h b/constraint_analysis/include/constraint_analysis/ca_parser.h
index 7bc41e60daf9d1c7263a7ce8918b933920ae8b02..55574447781e92f7df7c65183b6efb6b190465c2 100644
--- a/constraint_analysis/include/constraint_analysis/ca_parser.h
+++ b/constraint_analysis/include/constraint_analysis/ca_parser.h
@@ -20,6 +20,9 @@
* This file is part of UNICADO.
*/
+#ifndef SRC_CONSTRAINT_ANALYSIS_PARSER_H_
+#define SRC_CONSTRAINT_ANALYSIS_PARSER_H_
+
#include <cmath>
#include <vector>
#include <numbers>
@@ -115,4 +118,5 @@ public:
auto get_beta(const std::string segment_from) -> const double;
};
-//#endif // CONSTRAINTANALYSIS_CONSTRAINTANALYSIS_H_
\ No newline at end of file
+
+#endif
\ No newline at end of file
diff --git a/constraint_analysis/include/constraint_analysis/ca_plotting.h b/constraint_analysis/include/constraint_analysis/ca_plotting.h
index ff29d0f918ddcedfd8d292d4307ff0310a5c1b92..aaf5a3ffcf83e6491a41a8bd9bfb800780cfba80 100644
--- a/constraint_analysis/include/constraint_analysis/ca_plotting.h
+++ b/constraint_analysis/include/constraint_analysis/ca_plotting.h
@@ -20,6 +20,10 @@
* This file is part of UNICADO.
*/
+
+#ifndef SRC_CONSTRAINT_ANALYSIS_PLOTTING_H_
+#define SRC_CONSTRAINT_ANALYSIS_PLOTTING_H_
+
#include <cmath>
#include <vector>
#include <numbers>
@@ -111,4 +115,6 @@ public:
axis->area(x, y)->face_color(rgba_color);
legend.push_back("Infeasible Area");
}
-};
\ No newline at end of file
+};
+
+#endif
\ No newline at end of file
diff --git a/constraint_analysis/src/ca_functions.cpp b/constraint_analysis/src/ca_functions.cpp
index f297141bea91c64adc642c30870081269189d4c3..0262ec1f7f139ae4eb3f53824ec6b739cf1351d8 100644
--- a/constraint_analysis/src/ca_functions.cpp
+++ b/constraint_analysis/src/ca_functions.cpp
@@ -160,18 +160,22 @@ namespace Mattingly {
return T_over_W_data;
}
- auto constraint_analysis::takeoff_climb_angle(const std::vector<double>& W_over_S_data, double K_1, double K_2, double C_D0, double weight_fraction, double thrust_lapse, double Mach, double altitude, double CL_max_TO, double takeoff_climb_angle) -> const std::vector<double>{
+ auto constraint_analysis::takeoff_climb_angle(const std::vector<double>& W_over_S_data, const std::vector<double>& C_D_data, double weight_fraction, double thrust_lapse, double Mach, double altitude, double CL_max_TO, double takeoff_climb_angle) -> const std::vector<double>{
std::vector<double> T_over_W_data;
// Lambda function for calculating T_over_W
- auto calculate_T_over_W = [this](double W_over_S, double K_1, double K_2, double C_D0, double weight_fraction, double thrust_lapse, double Mach, double altitude, double CL_max_TO, double takeoff_climb_angle) {
- return (weight_fraction / thrust_lapse) * (K_1 * CL_max_TO / pow(k_T0, 2.0) + K_2 + (C_D0 / (CL_max_TO / pow(k_T0, 2.0))) + std::sin(takeoff_climb_angle*PI/180.0));
+ auto calculate_T_over_W = [this](double W_over_S, double C_D, double weight_fraction, double thrust_lapse, double Mach, double altitude, double CL_max_TO, double takeoff_climb_angle) {
+ return (weight_fraction / thrust_lapse) * (C_D/(CL_max_TO / pow(k_T0, 2.0)) + std::sin(takeoff_climb_angle * PI / 180.0));
};
// For each W_over_S calculate the T_over_W
- for (double W_over_S : W_over_S_data) {
- T_over_W_data.push_back(calculate_T_over_W(W_over_S, K_1, K_2, C_D0, weight_fraction, thrust_lapse, Mach, altitude, CL_max_TO, takeoff_climb_angle)); // Speichere Ergebnis in T_over_W_data
+ for (size_t i = 0; i < W_over_S_data.size(); ++i) {
+ double W_over_S = W_over_S_data[i];
+
+ double C_D = C_D_data[i];
+
+ T_over_W_data.push_back(calculate_T_over_W(W_over_S, C_D, weight_fraction, thrust_lapse, Mach, altitude, CL_max_TO, takeoff_climb_angle)); // Speichere Ergebnis in T_over_W_data
}
return T_over_W_data;
diff --git a/constraint_analysis/src/energy_based/energy_based.cpp b/constraint_analysis/src/energy_based/energy_based.cpp
index fc42ce6a9cf6e7b55ee1e4e577bcb4dc5c91cce3..b47ebd552647b3d6292cfaf4bfe76ec83fc75ba0 100644
--- a/constraint_analysis/src/energy_based/energy_based.cpp
+++ b/constraint_analysis/src/energy_based/energy_based.cpp
@@ -57,11 +57,23 @@ void EnergyBased::initialize()
void EnergyBased::operator()()
{
- /*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
- const int WS_start = 50;
- const int WS_end = 820;
- constexpr size_t WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
- std::vector<int> numbers(WS_size);
+ const std::string mode_selection(this->configuration()->at("module_configuration_file/program_settings/mode_selection/value"));
+
+ int WS_start = 0;
+ int WS_end = 0;
+ size_t WS_size = 1;
+ if (mode_selection == "mode_0")
+ {
+ /*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
+ WS_start = 50;
+ WS_end = 820;
+ WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
+ }
+ else if (mode_selection == "mode_1")
+ {
+ WS_start = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
+ WS_size = 1; // Constant W/S
+ }
std::vector<double> W_over_S_data(WS_size);
@@ -118,8 +130,8 @@ void EnergyBased::operator()()
/* Take-off Segment */
double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
double s_G_TO = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
- double takeoff_climb_angle = this->configuration()->at("module_configuration_file/program_settings/takeoff_climb_angle/value");
- double M_TO = this->configuration()->at("module_configuration_file/program_settings/Mach_TO/value");
+ double takeoff_climb_angle = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/takeoff_climb/takeoff_climb_angle/value");
+ double M_TO = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/takeoff_ground_roll/Mach_TO/value");
double thrust_lapse_TO = calculate_thrust_lapse(M_TO, 0.0);
double weight_fraction_TO = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
const auto to_ground_roll = constraint_analysis_tool.takeoff_ground_roll(
@@ -152,7 +164,7 @@ void EnergyBased::operator()()
/* OEI Climb */
double climb_speed = this->aircraft_xml()->at(
"aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/climb_speed_schedule/climb_speed_below_FL100/value");
- double climb_gradient_OEI = this->configuration()->at("module_configuration_file/program_settings/climb_gradient_OEI/value");
+ double climb_gradient_OEI = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/OEI_climb/climb_gradient_OEI/value");
double M_climb = climb_speed / atm.getSpeedOfSound(0.0);
double num_engines = get_num_engines();
@@ -229,7 +241,7 @@ void EnergyBased::operator()()
throw std::runtime_error("Invalid method for calculating CD");
}
- double minimum_climb_rate = this->configuration()->at("module_configuration_file/program_settings/minimum_climb_rate/value");
+ double minimum_climb_rate = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/SEP/minimum_climb_rate/value");
/* Service ceiling with the minimum climb rate requirement */
const auto service_ceiling = constraint_analysis_tool.constant_speed_climb(
@@ -258,8 +270,8 @@ void EnergyBased::operator()()
double AoA_2 = landing_config.interpalpha(CL_max_L - CL_max_L * 0.9);
double CL_alpha = (CL_max_L - (CL_max_L - CL_max_L * 0.9)) / ((AoA_1 - AoA_2) * PI / 180.0);
double V_TO_L = M_TO * atm.getSpeedOfSound(0.0);
- double w_g = this->configuration()->at("module_configuration_file/program_settings/gust_speed/value");
- double dn_G = this->configuration()->at("module_configuration_file/program_settings/gust_load_factor/value");
+ double w_g = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_speed/value");
+ double dn_G = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_load_factor/value");
const auto gust = constraint_analysis_tool.gust(CL_alpha, 0.0, V_TO_L, w_g, dn_G, weight_fraction_TO);
/*Define Constraint out of each non-vertical constraint case*/
@@ -281,8 +293,14 @@ void EnergyBased::operator()()
double safety_factor = this->configuration()->at("module_configuration_file/program_settings/safety_factor/value");
- /* Find design point with safety factor*/
- min_finder.find_design_point(safety_factor);
+ const double buffer_factor(this->configuration()->at("module_configuration_file/program_settings/buffer_factor/value"));
+
+ /*Get the previous design point*/
+ double previous_WS = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
+ double previous_TW = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/thrust_to_weight/value");
+
+ /*Find design point with safety factor*/
+ min_finder.find_converged_design_point(previous_TW, previous_WS * 9.81, safety_factor, buffer_factor);
/*Update design point in ACxml*/
min_finder.update_design_point();
diff --git a/constraint_analysis/src/energy_based/energy_based.h b/constraint_analysis/src/energy_based/energy_based.h
index 3b51508398af75626d834156c8bb815ea7da472e..d027e566ba063c9872e53f39ae5fecfbcc380293 100644
--- a/constraint_analysis/src/energy_based/energy_based.h
+++ b/constraint_analysis/src/energy_based/energy_based.h
@@ -20,8 +20,11 @@
* This file is part of UNICADO.
*/
+
+#ifndef SRC_CONSTRAINT_ANALYSIS_ENERGY_BASED_H_
+#define SRC_CONSTRAINT_ANALYSIS_ENERGY_BASED_H_
+
/* === Includes === */
-#include "../ca_strategy.h"
#include <filesystem>
#include <memory>
#include <optional>
@@ -29,9 +32,15 @@
#include <unordered_set>
#include <string>
#include <string_view>
+#include "../ca_strategy.h"
+#include <engine/engine.h>
+#include "../../include/constraint_analysis/ca_functions.h"
+#include "../../include/constraint_analysis/ca_parser.h"
+#include "../constraint_analysis.h"
+#include "../../include/constraint_analysis/ca_plotting.h"
+#include "../io/aircraft_xml.h"
+#include "../../include/constraint_analysis/ca_minfinder.h"
-#ifndef SRC_CONSTRAINT_ANALYSIS_ENERGY_BASED_H_
-#define SRC_CONSTRAINT_ANALYSIS_ENERGY_BASED_H_
class EnergyBased : public ConstraintAnalysisStrategy
{
@@ -58,6 +67,10 @@ public:
void operator()(); // NOLINT runtime/references
+ void assess_constraints(std::vector<double>& W_over_S_data, Mattingly::constraint_analysis& constraint_analysis_tool, readMission& mission_data, Engine& engine);
+
+ std::vector<Constraint> constraint_list = {};
+ std::vector<double> boundaries = {};
}; // namespace design
#endif // SRC_CONSTRAINT_ANALYSIS_ENERGY_BASED_H_
\ No newline at end of file
diff --git a/constraint_analysis/src/energy_based/energy_based_refactored.cpp b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb07a1cd11730fb3e131f7713c2684e6d2ecf249
--- /dev/null
+++ b/constraint_analysis/src/energy_based/energy_based_refactored.cpp
@@ -0,0 +1,347 @@
+/*
+ * UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+ *
+ * Copyright (C) 2024 UNICADO consortium
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <https://www.gnu.org/licenses/>.
+ *
+ * Description:
+ * This file is part of UNICADO.
+ */
+
+#include "energy_based.h"
+#include <algorithm>
+#include <filesystem>
+#include <format>
+#include <list>
+#include <ranges>
+#include <cmath>
+#include <vector>
+#include <numbers>
+#include <iostream>
+#include <string>
+#include <moduleBasics/module.h>
+#include <engine/engine.h>
+#include "../../include/constraint_analysis/ca_functions.h"
+#include "../../include/constraint_analysis/ca_parser.h"
+#include "../constraint_analysis.h"
+#include "../../include/constraint_analysis/ca_plotting.h"
+#include "../io/aircraft_xml.h"
+#include "../../include/constraint_analysis/ca_minfinder.h"
+
+
+void EnergyBased::initialize()
+{
+ if (!this->aircraft_xml()) {
+ throw std::runtime_error("[EnergyBased::initialize()] aircraft_xml is null!");
+ };
+ if (!this->configuration()) {
+ throw std::runtime_error("[EnergyBased::initialize()] config_xml is null!");
+ };
+ if (!this->polar()) {
+ throw std::runtime_error("[EnergyBased::initialize()] polar_xml is null!");
+ };
+}
+
+void EnergyBased::operator()()
+{
+
+ const std::string mode_selection(this->configuration()->at("module_configuration_file/program_settings/mode_selection/value"));
+
+ int WS_start = 0;
+ int WS_end = 0;
+ size_t WS_size = 1;
+
+ if (mode_selection == "mode_0")
+ {
+ /*Define relevant Wing_loading range. Selected values should be valid for all possible configurations*/
+ WS_start = 50;
+ WS_end = 820;
+ WS_size = WS_end - WS_start + 1; // Range from 50 to 820, inclusive
+ }
+ else if (mode_selection == "mode_1")
+ {
+ WS_start = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
+ WS_size = 1; // Constant W/S
+ }
+
+ std::vector<double> W_over_S_data(WS_size);
+
+ /*Transform from kg/m^2 into N/m^2*/
+ for (int i = 0; i < WS_size; ++i) {
+ W_over_S_data[i] = (WS_start + i) * 9.81;
+ }
+
+ Engine engine = Engine(this->engine_directory());
+
+ Mattingly::constraint_analysis constraint_analysis_tool;
+
+ atmosphere atm;
+
+ // Set atmosphere according to the delta ISA requirement for the sizing process
+ double temperature_ISA = ISA_TEMPERATURE + this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/delta_ISA/value");
+ atm.setAtmosphere(0.0, temperature_ISA, ISA_PRESSURE);
+
+ constraint_analysis_tool.atm = atm;
+
+ std::filesystem::path missionCSV = this->mission_directory();
+
+ for (const auto& entry : fs::directory_iterator(this->mission_directory())) {
+ std::string filename = entry.path().filename().string();
+
+ if (filename.find("design_mission") != std::string::npos && filename.find("_out.csv") != std::string::npos) {
+ missionCSV /= filename;
+ break;
+ }
+ }
+
+ readMission mission_data(missionCSV);
+
+ assess_constraints(W_over_S_data, constraint_analysis_tool, mission_data, engine);
+
+ /* Define Simple_Analysis*/
+ Simple_Analysis min_finder = Simple_Analysis(this->constraint_list, this->aircraft_xml(), this->boundaries);
+
+ /* Find dominant curve of the design chart*/
+ min_finder.find_dominant_curve();
+
+ double safety_factor = this->configuration()->at("module_configuration_file/program_settings/safety_factor/value");
+
+ const double buffer_factor(this->configuration()->at("module_configuration_file/program_settings/buffer_factor/value"));
+
+ /*Get the previous design point*/
+ double previous_WS = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/wing_loading/value");
+ double previous_TW = this->aircraft_xml()->at("aircraft_exchange_file/sizing_point/thrust_to_weight/value");
+
+ /*Find design point with safety factor*/
+ min_finder.find_converged_design_point(previous_TW, previous_WS * 9.81, safety_factor, buffer_factor);
+
+ /*Update design point in ACxml*/
+ min_finder.update_design_point();
+
+ std::string plot_dir_string = this->plot_dir().string();
+
+ // Add output of design point in logfile
+
+ myRuntimeInfo->out << "New wing loading: " << min_finder.design_wing_loading * 9.81 << " N/m^2. New thrust to weight ratio: " << min_finder.design_thrust_to_weight << std::endl;
+
+ min_finder.write_design_point(plot_dir_string + "/csv_files/design_point.csv");
+
+ /* Plot the complete design chart with constraints, boundaries and design point*/
+ std::string plotting = this->configuration()->at("module_configuration_file/control_settings/plot_output/enable/value");
+ if (plotting == "true") {
+ ConstraintAnalysis_Plot ca_plot;
+
+ ca_plot.fill_infeasible_area(W_over_S_data, min_finder.dominant_thrust_to_weight, boundaries);
+ for (auto ca_case : this->constraint_list)
+ {
+ ca_plot.add_curve(ca_case.wing_loading, ca_case.thrust_to_weight, ca_case.constraint_case);
+ }
+ for (auto boundary : this->boundaries)
+ {
+ ca_plot.add_curve(boundary, "W/S limit");
+ }
+ ca_plot.add_curve({ min_finder.design_wing_loading * 9.81 }, { min_finder.design_thrust_to_weight }, "Design Point");
+
+ ca_plot.save_figure(plot_dir_string + "/constraint_plot.svg");
+ };
+}
+
+void EnergyBased::assess_constraints(std::vector<double>& W_over_S_data, Mattingly::constraint_analysis& constraint_analysis_tool, readMission& mission_data, Engine& engine)
+{
+ auto calculate_thrust_lapse = [&constraint_analysis_tool, &engine](double M, double altitude) {
+ const double thrust_lapse = engine.get_thrust_lapse("cruise", constraint_analysis_tool.atm, M, altitude);
+ return thrust_lapse;
+ };
+
+ auto get_num_engines = [&]() {
+ double num_engines = 0;
+
+ // Get the node vector for all defined propulsion systems
+ auto propulsors = this->aircraft_xml()->getVector(
+ "aircraft_exchange_file/requirements_and_specifications/design_specification/propulsion/propulsor");
+
+ // Increment num_engines for each propulsor
+ for (const auto& propulsor : propulsors) {
+ num_engines += 1;
+ }
+ return num_engines;
+ };
+
+ auto std_constraints = this->configuration()->at(
+ "module_configuration_file/program_settings/constraint_selection/standard_set").getChildren();
+
+ for (auto constraint : std_constraints)
+ {
+ std::string activated = constraint->at("active/value");
+ if (activated == "true")
+ {
+ double CD = 0.0;
+ double CL = 0.0;
+ std::vector<double> CD_vector = {};
+ std::vector<double> CL_vector = {};
+
+ if (constraint->name == "takeoff_ground_roll")
+ {
+ double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
+ double s_G_TO = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
+ double M = constraint->at("Mach_TO/value");
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
+ const auto constraint_vector = constraint_analysis_tool.takeoff_ground_roll(
+ W_over_S_data,
+ weight_fraction,
+ thrust_lapse,
+ M,
+ 0.0,
+ s_G_TO,
+ CL_max_TO);
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "takeoff_climb")
+ {
+ double M = constraint->at("Mach_TO/value");
+ double takeoff_climb_angle = constraint->at("takeoff_climb_angle/value");
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
+ double CL_max_TO = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxT-O/value");
+ readPolar takeoff_config(this->polar(), "takeoff", "linear", M);
+ CD = takeoff_config.interpCD(CL_max_TO);
+ for (int i = 0; i < W_over_S_data.size(); ++i) {
+ CD_vector.push_back(CD);
+ };
+ const auto constraint_vector = constraint_analysis_tool.takeoff_climb_angle(
+ W_over_S_data,
+ CD_vector,
+ weight_fraction,
+ thrust_lapse,
+ M,
+ 0.0,
+ CL_max_TO,
+ takeoff_climb_angle);
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "OEI_climb")
+ {
+ /* OEI Climb */
+ double climb_speed = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/climb_speed_schedule/climb_speed_below_FL100/value");
+ double climb_gradient_OEI = constraint->at("climb_gradient_OEI/value");
+ double M = climb_speed / constraint_analysis_tool.atm.getSpeedOfSound(0.0);
+ double weight_fraction = mission_data.get_beta(" takeoff"); //from mission analysis
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double num_engines = get_num_engines();
+
+ readPolar climb_config(this->polar(), "takeoff", "linear", M);
+ for (int i = 0; i < W_over_S_data.size(); ++i) {
+ CD = climb_config.getCD(W_over_S_data[i], 1.0, weight_fraction, (climb_speed / constraint_analysis_tool.atm.getSpeedOfSound(0.0)), 0.0);
+ CD_vector.push_back(CD);
+ };
+
+ double thrust_factor_OEI = thrust_lapse * ((num_engines - 1) / num_engines);
+ double climb_rate_OEI = climb_gradient_OEI * climb_speed;
+ /*minimal climb rate of 2.4% defined in CS25 for second segment climb*/
+ const auto constraint_vector = constraint_analysis_tool.constant_speed_climb(
+ W_over_S_data,
+ CD_vector,
+ weight_fraction,
+ thrust_factor_OEI,
+ M,
+ 0.0,
+ climb_rate_OEI); /* DEFINE CLIMB RATE CONSTRAINT */
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "top_of_climb")
+ {
+
+ }
+ else if (constraint->name == "cruise")
+ {
+
+ }
+ else if (constraint->name == "SEP")
+ {
+ /* Service Ceiling */
+ double altitude = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/design_mission/initial_cruise_altitude/value");
+ double M = this->aircraft_xml()->at(
+ "aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/study_mission/initial_cruise_mach_number/value");
+
+ double weight_fraction = mission_data.get_beta(" climb_to_cruise");
+ readPolar clean_config_M_max(this->polar(), "clean", "linear", M);
+ double thrust_lapse = calculate_thrust_lapse(M, altitude);
+ CD_vector = {};
+ CL_vector = {};
+
+ for (int i = 0; i < W_over_S_data.size(); ++i) {
+ CD = clean_config_M_max.getCD(W_over_S_data[i], 1.0, weight_fraction, M, altitude);
+ CD_vector.push_back(CD);
+ }
+
+ double minimum_climb_rate = constraint->at("minimum_climb_rate/value");
+
+ /* Service ceiling with the minimum climb rate requirement */
+ const auto constraint_vector = constraint_analysis_tool.constant_speed_climb(
+ W_over_S_data,
+ CD_vector,
+ weight_fraction,
+ thrust_lapse,
+ M,
+ altitude,
+ minimum_climb_rate);
+ Constraint Case(W_over_S_data, constraint_vector, constraint->name);
+ this->constraint_list.push_back(Case);
+ }
+ else if (constraint->name == "approach_speed")
+ {
+
+ }
+ else if (constraint->name == "landing_field_length")
+ {
+ /* Landing */
+ double M = constraint->at("Mach/value");
+ double thrust_lapse = calculate_thrust_lapse(M, 0.0);
+ double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
+ double my_B = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/additional_requirements/landing_gear/braking_coefficient/value");
+ double s_G_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/takeoff_distance/value");
+ readPolar landing_config(this->polar(), "landing", "linear", M);
+ double CD_max_L = landing_config.interpCD(CL_max_L);
+ double max_landing_weight = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_landing_mass/mass_properties/mass/value");
+ double MTOM = this->aircraft_xml()->at("aircraft_exchange_file/analysis/masses_cg_inertia/maximum_takeoff_mass/mass_properties/mass/value");
+ double weight_fraction = max_landing_weight / MTOM;
+ const auto constraint_border = constraint_analysis_tool.braking_roll(CD_max_L, CL_max_L, weight_fraction, thrust_lapse, M, 0.0, my_B, s_G_L);
+ this->boundaries.push_back(constraint_border);
+ }
+ else if (constraint->name == "gust")
+ {
+ /* Gust */
+ double M = constraint->at("Mach/value");
+ readPolar landing_config(this->polar(), "landing", "linear", M);
+ double weight_fraction = mission_data.get_beta(" takeoff", 1.0); //from mission analysis
+ double CL_max_L = this->aircraft_xml()->at("aircraft_exchange_file/analysis/aerodynamics/lift_coefficients/C_LmaxLanding/value");
+ double AoA_1 = landing_config.interpalpha(CL_max_L);
+ double AoA_2 = landing_config.interpalpha(CL_max_L - CL_max_L * 0.9);
+ double CL_alpha = (CL_max_L - (CL_max_L - CL_max_L * 0.9)) / ((AoA_1 - AoA_2) * PI / 180.0);
+ double V_TO_L = this->aircraft_xml()->at("aircraft_exchange_file/requirements_and_specifications/requirements/top_level_aircraft_requirements/flight_envelope/maximum_approach_speed/value");
+ double w_g = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_speed/value");
+ double dn_G = this->configuration()->at("module_configuration_file/program_settings/constraint_selection/standard_set/gust/gust_load_factor/value");
+ const auto constraint_border = constraint_analysis_tool.gust(CL_alpha, 0.0, V_TO_L, w_g, dn_G, weight_fraction);
+ this->boundaries.push_back(constraint_border);
+ }
+ }
+ };
+};
\ No newline at end of file