From b4b28f3407d466207b5c30840459e3a00caa070d Mon Sep 17 00:00:00 2001 From: Jasper Bussemaker Date: Wed, 18 Feb 2026 08:23:21 +0100 Subject: [PATCH 1/6] Replace np.row_stack with np.vstack --- sb_arch_opt/algo/arch_sbo/infill.py | 2 +- sb_arch_opt/algo/egor_interface/algo.py | 6 +++--- sb_arch_opt/algo/segomoe_interface/algo.py | 14 +++++++------- sb_arch_opt/algo/tpe_interface/api.py | 2 +- sb_arch_opt/design_space.py | 4 ++-- sb_arch_opt/pareto_front.py | 4 ++-- sb_arch_opt/problems/gnc.py | 6 +++--- sb_arch_opt/problems/rocket.py | 2 +- sb_arch_opt/sampling.py | 14 +++++++------- 9 files changed, 27 insertions(+), 27 deletions(-) diff --git a/sb_arch_opt/algo/arch_sbo/infill.py b/sb_arch_opt/algo/arch_sbo/infill.py index dee71ef..2f8c4ec 100644 --- a/sb_arch_opt/algo/arch_sbo/infill.py +++ b/sb_arch_opt/algo/arch_sbo/infill.py @@ -283,7 +283,7 @@ def _g(_): else: x_optimized.append(x_ref_i) - return Population.new(X=np.row_stack(x_optimized)) + return Population.new(X=np.vstack(x_optimized)) @staticmethod def get_pareto_front(f: np.ndarray) -> np.ndarray: diff --git a/sb_arch_opt/algo/egor_interface/algo.py b/sb_arch_opt/algo/egor_interface/algo.py index 4543abb..bd829ef 100644 --- a/sb_arch_opt/algo/egor_interface/algo.py +++ b/sb_arch_opt/algo/egor_interface/algo.py @@ -217,9 +217,9 @@ def _run_infills(self, n_infills: int): x, x_failed, y = self._get_xy(pop) # Update - self._x = np.row_stack([self._x, x]) - self._y = np.row_stack([self._y, y]) - self._x_failed = np.row_stack([self._x_failed, x_failed]) + self._x = np.vstack([self._x, x]) + self._y = np.vstack([self._y, y]) + self._x_failed = np.vstack([self._x_failed, x_failed]) # Store results if self._results_folder is not None: diff --git a/sb_arch_opt/algo/segomoe_interface/algo.py b/sb_arch_opt/algo/segomoe_interface/algo.py index 9354982..76a7cda 100644 --- a/sb_arch_opt/algo/segomoe_interface/algo.py +++ b/sb_arch_opt/algo/segomoe_interface/algo.py @@ -316,9 +316,9 @@ def _dummy_f_grouped(_): return sego.get_x(i=-1) def _tell_infill(self, x, x_failed, y): - self._x = np.row_stack([self._x, x]) if self._x is not None else x - self._y = np.row_stack([self._y, y]) if self._y is not None else y - self._x_failed = np.row_stack([self._x_failed, x_failed]) if self._x_failed is not None else x_failed + self._x = np.vstack([self._x, x]) if self._x is not None else x + self._y = np.vstack([self._y, y]) if self._y is not None else y + self._x_failed = np.vstack([self._x_failed, x_failed]) if self._x_failed is not None else x_failed self._save_results() def _get_sego(self, f_grouped): @@ -465,10 +465,10 @@ def get_population(self, x: np.ndarray, y: np.ndarray, x_failed: np.ndarray = No f, g, h = self._split_y(y) if x_failed is not None and len(x_failed) > 0: - x = np.row_stack([x, x_failed]) - f = np.row_stack([f, np.zeros((x_failed.shape[0], f.shape[1]))*np.inf]) - g = np.row_stack([g, np.zeros((x_failed.shape[0], g.shape[1]))*np.inf]) - h = np.row_stack([h, np.zeros((x_failed.shape[0], h.shape[1]))*np.inf]) + x = np.vstack([x, x_failed]) + f = np.vstack([f, np.zeros((x_failed.shape[0], f.shape[1]))*np.inf]) + g = np.vstack([g, np.zeros((x_failed.shape[0], g.shape[1]))*np.inf]) + h = np.vstack([h, np.zeros((x_failed.shape[0], h.shape[1]))*np.inf]) kwargs = {'X': x, 'F': f, 'G': g, 'H': h} pop = Population.new(**kwargs) diff --git a/sb_arch_opt/algo/tpe_interface/api.py b/sb_arch_opt/algo/tpe_interface/api.py index a6cea7c..6325f33 100644 --- a/sb_arch_opt/algo/tpe_interface/api.py +++ b/sb_arch_opt/algo/tpe_interface/api.py @@ -170,7 +170,7 @@ def __init__(self): super().__init__(sampling=None) def do(self, problem, n_samples, **kwargs): - x_init = np.row_stack([self.interface.ask_init() for _ in range(n_samples)]) + x_init = np.vstack([self.interface.ask_init() for _ in range(n_samples)]) return Population.new(X=x_init) diff --git a/sb_arch_opt/design_space.py b/sb_arch_opt/design_space.py index 7f8056c..7e1ca2e 100644 --- a/sb_arch_opt/design_space.py +++ b/sb_arch_opt/design_space.py @@ -647,8 +647,8 @@ def _get_all_discrete_x_by_trial_and_imputation(self): x_repair = x_repair[is_not_repaired, :] is_active = is_active[is_not_repaired, :] - x_discr = np.row_stack([x_discr, x_repair]) - is_act_discr = np.row_stack([is_act_discr, is_active.astype(bool)]) + x_discr = np.vstack([x_discr, x_repair]) + is_act_discr = np.vstack([is_act_discr, is_active.astype(bool)]) # Impute continuous values self.impute_x(x_discr, is_act_discr) diff --git a/sb_arch_opt/pareto_front.py b/sb_arch_opt/pareto_front.py index db78596..9751443 100644 --- a/sb_arch_opt/pareto_front.py +++ b/sb_arch_opt/pareto_front.py @@ -140,9 +140,9 @@ def _calc_pareto_set_front(self, *_, pop_size=None, n_gen_min=10, n_repeat=4, n_ ps = res.X pf = res.F else: - pf_merged = np.row_stack([pf, res.F]) + pf_merged = np.vstack([pf, res.F]) i_non_dom = NonDominatedSorting().do(pf_merged, only_non_dominated_front=True) - ps = np.row_stack([ps, res.X])[i_non_dom, :] + ps = np.vstack([ps, res.X])[i_non_dom, :] pf = pf_merged[i_non_dom, :] # Reduce size of Pareto front to a predetermined amount to ease Pareto-front-related calculations diff --git a/sb_arch_opt/problems/gnc.py b/sb_arch_opt/problems/gnc.py index 6f4ef53..c3f11c1 100644 --- a/sb_arch_opt/problems/gnc.py +++ b/sb_arch_opt/problems/gnc.py @@ -175,8 +175,8 @@ def _gen_all_discrete_x(self) -> Optional[Tuple[np.ndarray, np.ndarray]]: x_rows.append(x_combs) is_active_rows.append(is_act_combs) - x_all = np.row_stack(x_rows) - is_active_all = np.row_stack(is_active_rows) + x_all = np.vstack(x_rows) + is_active_all = np.vstack(is_active_rows) return x_all, is_active_all def _get_discrete_x_combs_type(self, x_base, j, n_objs): @@ -234,7 +234,7 @@ def _iter_conns(n_src_, n_tgt_): n_combinations += 1 if return_conns: - n_combinations = np.row_stack(n_combinations) + n_combinations = np.vstack(n_combinations) n_comb_conn[n_src, n_tgt] = n_comb_conn[n_tgt, n_src] = n_combinations if return_conns: diff --git a/sb_arch_opt/problems/rocket.py b/sb_arch_opt/problems/rocket.py index 6656e91..61967dd 100644 --- a/sb_arch_opt/problems/rocket.py +++ b/sb_arch_opt/problems/rocket.py @@ -166,7 +166,7 @@ def _gen_all_discrete_x(self) -> Optional[Tuple[np.ndarray, np.ndarray]]: x_stage[:, 1:1+x_engines.shape[1]] = x_engines x_stages.append(x_stage) - x_stages = np.row_stack(x_stages) + x_stages = np.vstack(x_stages) x_all = np.repeat(x_stages, 3, axis=0) x_all[:, [11]] = np.tile(np.array([np.arange(3)]).T, (x_stages.shape[0], 1)) # Head shape diff --git a/sb_arch_opt/sampling.py b/sb_arch_opt/sampling.py index de37507..25ff845 100644 --- a/sb_arch_opt/sampling.py +++ b/sb_arch_opt/sampling.py @@ -339,8 +339,8 @@ def _choice(n_choose, n_from, replace=True): is_active.append(is_act_all[x_all_choose, :]) i_x_sampled[x_all_choose] = True - x = np.row_stack(x) - is_active = np.row_stack(is_active) + x = np.vstack(x) + is_active = np.vstack(is_active) # Uniformly add discrete vectors if there are not enough (can happen if some groups are very small and there # are no continuous dimensions) @@ -354,8 +354,8 @@ def _choice(n_choose, n_from, replace=True): else: i_from_group = np.arange(x_available.shape[0]) - x = np.row_stack([x, x_available[i_from_group, :]]) - is_active = np.row_stack([is_active, is_act_available[i_from_group, :]]) + x = np.vstack([x, x_available[i_from_group, :]]) + is_active = np.vstack([is_active, is_act_available[i_from_group, :]]) return x, is_active @@ -489,7 +489,7 @@ def _choice(n_choose, n_from, replace=True): i_opt_sampled = _choice(n_add, len(opt_values[i_dv])) x_add[:, i_dv] = opt_values[i_dv][i_opt_sampled] - x = x_add if x is None else np.row_stack([x, x_add]) + x = x_add if x is None else np.vstack([x, x_add]) # Correct and remove duplicates x, is_active = self._correct(problem, repair, x) @@ -509,8 +509,8 @@ def _choice(n_choose, n_from, replace=True): if x.shape[0] < n_samples and has_x_cont: n_add = n_samples-x.shape[0] i_select_dup = _choice(n_add, x.shape[0]) - x = np.row_stack(x, x[i_select_dup, :]) - is_active = np.row_stack(is_active, is_active[i_select_dup, :]) + x = np.vstack(x, x[i_select_dup, :]) + is_active = np.vstack(is_active, is_active[i_select_dup, :]) return x, is_active From 168bb5d9a85a9feca2d3936481d451787a9592ca Mon Sep 17 00:00:00 2001 From: Jasper Bussemaker Date: Wed, 18 Feb 2026 10:31:29 +0100 Subject: [PATCH 2/6] Version bump --- sb_arch_opt/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sb_arch_opt/__init__.py b/sb_arch_opt/__init__.py index bb64aa4..4a9b978 100644 --- a/sb_arch_opt/__init__.py +++ b/sb_arch_opt/__init__.py @@ -1 +1 @@ -__version__ = '1.6.1' +__version__ = '1.6.2' From 853e200a50dc2e58a7d15e22598368128375fc72 Mon Sep 17 00:00:00 2001 From: Jasper Bussemaker Date: Wed, 22 Jul 2026 09:39:33 +0200 Subject: [PATCH 3/6] Version bump --- sb_arch_opt/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sb_arch_opt/__init__.py b/sb_arch_opt/__init__.py index 4a9b978..4574cc8 100644 --- a/sb_arch_opt/__init__.py +++ b/sb_arch_opt/__init__.py @@ -1 +1 @@ -__version__ = '1.6.2' +__version__ = '1.6.3' From bc17f90bba49f2de7a9ce8cbe3fdebe7e8806f36 Mon Sep 17 00:00:00 2001 From: Jasper Bussemaker Date: Wed, 22 Jul 2026 09:41:56 +0200 Subject: [PATCH 4/6] Add optional constraint normalization and failure if orbit cannot be reached to the rocket problem --- sb_arch_opt/problems/rocket.py | 16 ++++++++++---- sb_arch_opt/problems/rocket_eval.py | 33 +++++++++++++++++++++++------ 2 files changed, 38 insertions(+), 11 deletions(-) diff --git a/sb_arch_opt/problems/rocket.py b/sb_arch_opt/problems/rocket.py index 61967dd..8594688 100644 --- a/sb_arch_opt/problems/rocket.py +++ b/sb_arch_opt/problems/rocket.py @@ -51,6 +51,8 @@ class RocketArch(HierarchyProblemBase): _head_shapes = [HeadShape.CONE, HeadShape.ELLIPTICAL, HeadShape.SPHERE] _less_constrained = False + normalized_constraints = True + fail_if_cannot_reach_orbit = False def __init__(self): check_dependency() @@ -85,7 +87,13 @@ def _arch_evaluate(self, x: np.ndarray, is_active_out: np.ndarray, f_out: np.nda lc = self._less_constrained rockets = self._get_rockets(x) for i, rocket in enumerate(rockets): - perf = RocketEvaluator.evaluate(rocket) + + perf = RocketEvaluator.evaluate( + rocket, + normalize_constraints=self.normalized_constraints, + fail_if_cannot_reach_orbit=self.fail_if_cannot_reach_orbit, + ) + f_out[i, :] = (np.log10(perf.cost), -np.log10(max(1., perf.payload_mass))) if lc: g_out[i, :] = (perf.delta_structural, perf.delta_payload, perf.delta_delta_v) @@ -221,15 +229,15 @@ def __repr__(self): from pymoo.core.population import Population from sb_arch_opt.sampling import HierarchicalSampling - # problem = RocketArch() + problem = RocketArch() # x_pf = problem.pareto_set() # f_pf = problem.pareto_front() # problem = LCRocketArch() # problem = SOLCRocketArch(obj=RocketObj.OBJ_COST) - problem = SOLCRocketArch(obj=RocketObj.OBJ_PAYLOAD) + # problem = SOLCRocketArch(obj=RocketObj.OBJ_PAYLOAD) # problem = SOLCRocketArch(obj=RocketObj.OBJ_WEIGHTED) - problem.plot_pf() + # problem.plot_pf() # f_pf = problem.pareto_front() # f_so = f_pf[:, 0] + f_pf[:, 1] diff --git a/sb_arch_opt/problems/rocket_eval.py b/sb_arch_opt/problems/rocket_eval.py index 9c95568..b0f7538 100644 --- a/sb_arch_opt/problems/rocket_eval.py +++ b/sb_arch_opt/problems/rocket_eval.py @@ -93,7 +93,7 @@ class RocketEvaluator: _rho_interp = None @classmethod - def evaluate(cls, rocket: Rocket) -> RocketPerformance: + def evaluate(cls, rocket: Rocket, normalize_constraints=False, fail_if_cannot_reach_orbit=False) -> RocketPerformance: check_dependency() # Calculate geometrical data @@ -164,7 +164,7 @@ def evaluate(cls, rocket: Rocket) -> RocketPerformance: length_ratio = rocket.ellipse_l_ratio if rocket.head_shape == HeadShape.ELLIPTICAL else 0 payload_mass, h_vector, v_vector, delta_delta_v = cls.calculate_trajectory( cone_angle, length_ratio, diameter, stage_thrusts, stage_structural_masses, stage_prop_masses, stage_mdots, - rocket.orbit_altitude, m_payload_fix=rocket.payload_mass) + rocket.orbit_altitude, m_payload_fix=rocket.payload_mass, normalize_constraint=normalize_constraints) # Cost estimation cost = cls.calculate_cost(stage_n_engines, stage_engine_masses, stage_solid_prop_masses, stage_h2_masses, @@ -174,6 +174,17 @@ def evaluate(cls, rocket: Rocket) -> RocketPerformance: delta_structural = cls.calculate_max_q_constraint(rocket.max_q, h_vector, v_vector) delta_payload = cls.calculate_payload_constraint(volume_available, payload_mass, rocket.payload_density) + if normalize_constraints: + delta_structural /= rocket.max_q + delta_payload /= volume_available + + if fail_if_cannot_reach_orbit and delta_delta_v < 0: + cost = math.nan + payload_mass = math.nan + delta_structural = math.nan + delta_payload = math.nan + delta_delta_v = math.nan + return RocketPerformance( cost=cost, payload_mass=payload_mass, delta_structural=delta_structural, delta_payload=delta_payload, delta_delta_v=-delta_delta_v, @@ -390,12 +401,13 @@ def modified_atmosphere(cls, x): @classmethod def calculate_trajectory(cls, cone_angle, length_ratio, diameter, T_stages, m_structural_stages, mp_stages, - mdot_stages, h_orbit_target, m_payload_fix=None): + mdot_stages, h_orbit_target, m_payload_fix=None, normalize_constraint=False): """Calculation of the launcher trajectory.""" check_dependency() mu = 3.986004418e14 r_earth = 6378e3 + v_orbit = (mu / (r_earth + h_orbit_target)) ** 0.5 # Drag coefficient calculation depending on head shape if cone_angle > 0: @@ -577,27 +589,34 @@ def stage_state_eq(t_, y): # Trajectory equation return v_orbit_final, h_vector_, v_vector_ def try_payload(m_payload): - v_orbit = (mu / (r_earth + h_orbit_target)) ** 0.5 try: v_final_, h_vector_, v_vector_ = simulate_trajectory(m_payload) # Orbit minimum speed v_target_diff = v_final_ - v_orbit + if normalize_constraint: + v_target_diff /= v_orbit + return v_target_diff, v_final_, h_vector_, v_vector_ except (IndexError, ValueError): - return -v_orbit, 0, [], [] + + v_diff = -v_orbit + if normalize_constraint: + v_diff /= v_orbit + + return v_diff, 0, [], [] # Evaluate for fixed payload mass if m_payload_fix: v_tgt_diff, _, h_vector, v_vector = try_payload(m_payload_fix) if v_tgt_diff < 0: - return 0, [], [] + return 0, [], [], v_tgt_diff return m_payload_fix, h_vector, v_vector, v_tgt_diff # Check if rocket could be feasible even without payload v_tgt_diff, _, _, _ = try_payload(0) - if v_tgt_diff <= 0: + if v_tgt_diff < 0: return 0, [], [], v_tgt_diff m_payload, res = opt.newton( From aea6dae92cb898ebaaaab688e0b4448074dcb27d Mon Sep 17 00:00:00 2001 From: Jasper Bussemaker Date: Wed, 22 Jul 2026 10:32:53 +0200 Subject: [PATCH 5/6] Allow converging to negative payloads --- sb_arch_opt/problems/rocket_eval.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/sb_arch_opt/problems/rocket_eval.py b/sb_arch_opt/problems/rocket_eval.py index b0f7538..18c9eea 100644 --- a/sb_arch_opt/problems/rocket_eval.py +++ b/sb_arch_opt/problems/rocket_eval.py @@ -495,6 +495,10 @@ def stage_state_eq(t_, y): # Trajectory equation h_first += h.tolist() v_first += v.tolist() + if len(pos) == 0: + # No second stage available + return v_first[-1], h_first, v_first + # Second stage h_0 = h[pos[0]] v_0 = v[pos[0]] @@ -562,6 +566,10 @@ def stage_state_eq(t_, y): # Trajectory equation h_second += h.tolist() v_second += v.tolist() + if len(pos2) == 0: + # No third stage available + return v_second[-1], h_first+h_second, v_first+v_second + # Third stage v_0 = v[pos2[0]] t_0 = t[pos2[0]] @@ -616,8 +624,8 @@ def try_payload(m_payload): # Check if rocket could be feasible even without payload v_tgt_diff, _, _, _ = try_payload(0) - if v_tgt_diff < 0: - return 0, [], [], v_tgt_diff + # if v_tgt_diff < 0: + # return 0, [], [], v_tgt_diff m_payload, res = opt.newton( lambda mp_: try_payload(mp_)[0], 100, tol=1., maxiter=50, full_output=True, disp=False) From 5754b287ac72528e2141b2f5a9dac1e33bd6c87e Mon Sep 17 00:00:00 2001 From: Jasper Bussemaker Date: Mon, 17 Aug 2026 15:27:53 +0200 Subject: [PATCH 6/6] Update installation instructions in readme --- README.md | 4 ++-- docs/index.md | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 039f0a9..35113a2 100644 --- a/README.md +++ b/README.md @@ -37,13 +37,13 @@ The library provides: First, create a conda environment (skip if you already have one): ``` -conda create --name opt python=3.11 +conda create --name opt python=3.12 conda activate opt ``` Then install the package: ``` -conda install "numpy<2.0" +conda install numpy pip install sb-arch-opt ``` diff --git a/docs/index.md b/docs/index.md index 213700d..d101b7b 100644 --- a/docs/index.md +++ b/docs/index.md @@ -31,13 +31,13 @@ The library provides: First, create a conda environment (skip if you already have one): ``` -conda create --name opt python=3.11 +conda create --name opt python=3.12 conda activate opt ``` Then install the package: ``` -conda install "numpy<2.0" +conda install numpy pip install sb-arch-opt ```