From 23f53bc9155dd1789aa04d665085c17ec3b27640 Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 14:30:50 +0100 Subject: [PATCH 01/55] Updated requirements --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 04ea1c0..7c245b5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -pyaedt +pyaedt==0.11.2 pyvista[jupyter] trame matplotlib \ No newline at end of file From 27a7392ac11ec108d521444041dda3c9f45f683f Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 15:09:50 +0100 Subject: [PATCH 02/55] Class constructors --- design.py | 35 +++++++++++++++++++++++++---------- generate_designs.ipynb | 24 ++++++++++++++---------- 2 files changed, 39 insertions(+), 20 deletions(-) diff --git a/design.py b/design.py index d5aee45..51fb8f4 100644 --- a/design.py +++ b/design.py @@ -3,11 +3,33 @@ from ansys.aedt.core import Desktop, Maxwell2d class Design: - def __init__(self, m2d=None): + def __init__(self, m2d): self.set_parameters() - if m2d is not None: - self.m2d = m2d + self.m2d = m2d + + @classmethod + def create(cls, project_name, design_name, file_name, **kwargs): + m2d = Maxwell2d( + project=project_name, + design=design_name, + solution_type="TransientXY", + **kwargs + ) + obj = cls(m2d) + obj.create_stator() + obj.save_project(file_name) + return obj + + @classmethod + def load(cls, file_name, **kwargs): + desktop = Desktop(**kwargs) + desktop.load_project(file_name) + + m2d = Maxwell2d() + m2d.set_active_design("Design01") + return cls(m2d) + def set_parameters(self): #materials self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" @@ -75,13 +97,6 @@ def set_parameters(self): self.rotor_r_min = self.mm_to_str('geom_params', 'DiaShaft')/2 self.rotor_r_max = self.mm_to_str('geom_params', 'DiaStatorGap')/2 - self.mm_to_str('geom_params', 'Airgap') - def load_stator(self, file_name, **kwargs): - desktop = Desktop(**kwargs) - desktop.load_project(file_name) - m2d = Maxwell2d() - m2d.set_active_design("Design01") - self.m2d = m2d - def create_stator(self): m2d = self.m2d modeler = m2d.modeler diff --git a/generate_designs.ipynb b/generate_designs.ipynb index 4a352e7..6a2ca34 100644 --- a/generate_designs.ipynb +++ b/generate_designs.ipynb @@ -9,7 +9,6 @@ "import os\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from ansys.aedt.core import Maxwell2d\n", "from design import Design\n", "from geometry import create_barrier\n", "\n", @@ -31,16 +30,21 @@ "CLS_EXIT = True #close on exit\n", "\n", "if not os.path.exists(file_name_aedt):\n", - " m2d = Maxwell2d(project=project_name, design=design_name, version=AEDT_VERSION,\n", - " non_graphical=NG_MODE, new_desktop=False,\n", - " close_on_exit=CLS_EXIT, student_version=False, solution_type=\"TransientXY\")\n", - " design = Design(m2d)\n", - " design.create_stator()\n", - " design.save_project(file_name_aedt)\n", + " design = Design.create(\n", + " project_name, design_name, file_name_aedt,\n", + " version=AEDT_VERSION,\n", + " non_graphical=NG_MODE,\n", + " new_desktop=False,\n", + " close_on_exit=CLS_EXIT,\n", + " )\n", "else:\n", - " design = Design()\n", - " design.load_stator(file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE,\n", - " new_desktop=False, close_on_exit=CLS_EXIT, student_version=False)" + " design = Design.load(\n", + " file_name_aedt,\n", + " version=AEDT_VERSION,\n", + " non_graphical=NG_MODE,\n", + " new_desktop=False,\n", + " close_on_exit=CLS_EXIT,\n", + " )" ] }, { From 6417965c074e8d38991a549d862dfaea35c4babb Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 15:16:49 +0100 Subject: [PATCH 03/55] Added typing --- design.py | 42 ++++++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/design.py b/design.py index 51fb8f4..cd10871 100644 --- a/design.py +++ b/design.py @@ -1,14 +1,15 @@ import os import numpy as np from ansys.aedt.core import Desktop, Maxwell2d +from ansys.aedt.core.modeler.modeler_2d import Modeler2D class Design: - def __init__(self, m2d): + def __init__(self, m2d: Maxwell2d) -> None: self.set_parameters() self.m2d = m2d @classmethod - def create(cls, project_name, design_name, file_name, **kwargs): + def create(cls, project_name: str, design_name: str, file_name: str, **kwargs) -> "Design": m2d = Maxwell2d( project=project_name, design=design_name, @@ -21,7 +22,7 @@ def create(cls, project_name, design_name, file_name, **kwargs): return obj @classmethod - def load(cls, file_name, **kwargs): + def load(cls, file_name: str, **kwargs) -> "Design": desktop = Desktop(**kwargs) desktop.load_project(file_name) @@ -30,7 +31,7 @@ def load(cls, file_name, **kwargs): return cls(m2d) - def set_parameters(self): + def set_parameters(self) -> None: #materials self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" @@ -97,9 +98,10 @@ def set_parameters(self): self.rotor_r_min = self.mm_to_str('geom_params', 'DiaShaft')/2 self.rotor_r_max = self.mm_to_str('geom_params', 'DiaStatorGap')/2 - self.mm_to_str('geom_params', 'Airgap') - def create_stator(self): + def create_stator(self) -> None: m2d = self.m2d modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) # Define design variables from the created dictionaries. modeler.model_units = "mm" @@ -452,8 +454,11 @@ def create_stator(self): plot_name=v, ) - def add_rotor(self): - rotor_id = self.m2d.modeler.create_polyline( + def add_rotor(self) -> None: + modeler = self.m2d.modeler + assert isinstance(modeler, Modeler2D) + + rotor_id = modeler.create_polyline( points=self.rot_points, segment_type=[ "Arc","Line", "Arc"], cover_surface=True, name="Rotor" ) self.rotor_id = rotor_id @@ -461,9 +466,9 @@ def add_rotor(self): rotor_id.color = (192, 192, 192) # rgb rotor_id.transparency = 0.0 - def add_rotor_barriers(self, barrier_points, n=1, segment_type=None, tol=1e-10): - m2d = self.m2d - modeler = m2d.modeler + def add_rotor_barriers(self, barrier_points, n: int = 1, segment_type=None, tol: float = 1e-10) -> None: + modeler = self.m2d.modeler + assert isinstance(modeler, Modeler2D) barrier_points = np.array(barrier_points) if np.array_equal(barrier_points[0], barrier_points[-1]): @@ -493,7 +498,7 @@ def add_rotor_barriers(self, barrier_points, n=1, segment_type=None, tol=1e-10): self.rotor_id.subtract(barrier_id) modeler.delete(barrier_id) - def compute(self, NUM_CORES=1): + def compute(self, NUM_CORES: int = 1): m2d = self.m2d m2d.mesh.assign_length_mesh( @@ -514,24 +519,25 @@ def compute(self, NUM_CORES=1): ) return solutions.data_magnitude() - def delete_rotor(self): + def delete_rotor(self) -> None: + assert isinstance(self.m2d.modeler, Modeler2D) self.m2d.modeler.delete(self.rotor_id) - def save_design(self, file_name, **kwargs): + def save_design(self, file_name: str, **kwargs) -> None: show = kwargs.pop('show', False) view = kwargs.pop('view', 'xy') self.m2d.plot(show=show, output_file=file_name, view=view) - def save_project(self, file_name=None): + def save_project(self, file_name: str | None = None) -> None: if file_name is None: self.m2d.save_project() else: self.m2d.save_project(file_name) - def close_project(self): + def close_project(self) -> None: self.m2d.close_desktop() - def analyze_results(self, Tor): + def analyze_results(self, Tor: np.ndarray) -> tuple[float, float, float]: TorAvg = np.mean(Tor[:-1]) TorAvgAC = np.mean(np.abs(Tor[:-1]-TorAvg)) TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) @@ -540,13 +546,13 @@ def analyze_results(self, Tor): return TorAvg, TorRmsAC, TorRippleRms - def mm_to_str(self, var, field): + def mm_to_str(self, var, field) -> float: val = getattr(self, var)[field] if not val.endswith('mm'): raise Exception('val must end with mm') return float(val[:-2]) - def print_results(self, TorAvg, TorRmsAC, TorRippleRms): + def print_results(self, TorAvg: float, TorRmsAC: float, TorRippleRms: float) -> None: print("\nTorque mean value: {:.2f} Nm".format(TorAvg)) # print("\nTorque ripple mean value: {:.2f} Nm".format(TorAvgAC)) print("\nTorque ripple rms value: {:.2f} Nm".format(TorRmsAC)) From 218945f5e3c424d597fc1a893a957a1e058eca21 Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 15:31:20 +0100 Subject: [PATCH 04/55] Simplify add_rotor_barrier --- design.py | 41 ++++++++++++++++------------------------- 1 file changed, 16 insertions(+), 25 deletions(-) diff --git a/design.py b/design.py index cd10871..e588856 100644 --- a/design.py +++ b/design.py @@ -466,40 +466,31 @@ def add_rotor(self) -> None: rotor_id.color = (192, 192, 192) # rgb rotor_id.transparency = 0.0 - def add_rotor_barriers(self, barrier_points, n: int = 1, segment_type=None, tol: float = 1e-10) -> None: + def add_rotor_barrier(self, barrier_points, segment_type=None) -> None: modeler = self.m2d.modeler assert isinstance(modeler, Modeler2D) - barrier_points = np.array(barrier_points) - if np.array_equal(barrier_points[0], barrier_points[-1]): - if np.linalg.norm(barrier_points[0] - barrier_points[-1]) < tol: - mid = 0.5*(barrier_points[0]+barrier_points[-1]) - barrier_points[0] = mid - barrier_points[-1] = mid - else: - raise Exception('First and last points are not equal') - barrier_points = np.round(barrier_points, 6) + # Round it for Ansys + barrier_points = np.round(barrier_points, 6) + + # Potentially add the z axis if barrier_points.shape[1] == 2: barrier_points = np.hstack((barrier_points, np.zeros((len(barrier_points),1)))) + + # Convert them into a string format and interpolate + points_str = [[str(y) for y in x] for x in barrier_points] + barrier_id = modeler.create_polyline( + points=points_str, segment_type=segment_type, cover_surface=True, name="Barrier" + ) - for i in range(n): - r = np.linalg.norm(barrier_points, axis=1) - r_proportion = (r - self.rotor_r_min) / (self.rotor_r_max - self.rotor_r_min) - r_min_new = i/n - r_max_new = (i+1)/n - r_proportion_new = r_min_new + (r_max_new - r_min_new) * r_proportion - r_new = self.rotor_r_min + (self.rotor_r_max - self.rotor_r_min) * r_proportion_new - barrier_points_new = barrier_points / np.reshape(r / r_new, (-1, 1)) - - qwe = [[str(y) for y in x] for x in barrier_points_new] - barrier_id = modeler.create_polyline( - points=qwe, segment_type=segment_type, cover_surface=True, name="Barrier" - ) - self.rotor_id.subtract(barrier_id) - modeler.delete(barrier_id) + # Remove the barrier + self.rotor_id.subtract(barrier_id) + modeler.delete(barrier_id) def compute(self, NUM_CORES: int = 1): m2d = self.m2d + assert m2d.mesh is not None + assert m2d.post is not None m2d.mesh.assign_length_mesh( assignment=self.rotor_id, From 1c7344861851f6ac4bac2e1613d2d4865b64c9c6 Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 17:11:46 +0100 Subject: [PATCH 05/55] Refactor --- design.py | 16 ------- generate_designs.ipynb | 61 ++++++-------------------- generators.py | 98 ++++++++++++++++++++++++++++++++++++++++++ geometry.py | 67 +++++++++++------------------ 4 files changed, 138 insertions(+), 104 deletions(-) create mode 100644 generators.py diff --git a/design.py b/design.py index e588856..5913dc2 100644 --- a/design.py +++ b/design.py @@ -527,25 +527,9 @@ def save_project(self, file_name: str | None = None) -> None: def close_project(self) -> None: self.m2d.close_desktop() - - def analyze_results(self, Tor: np.ndarray) -> tuple[float, float, float]: - TorAvg = np.mean(Tor[:-1]) - TorAvgAC = np.mean(np.abs(Tor[:-1]-TorAvg)) - TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) - TorRippleAvg = TorAvgAC/TorAvg*100 - TorRippleRms = TorRmsAC/TorAvg*100 - - return TorAvg, TorRmsAC, TorRippleRms def mm_to_str(self, var, field) -> float: val = getattr(self, var)[field] if not val.endswith('mm'): raise Exception('val must end with mm') return float(val[:-2]) - - def print_results(self, TorAvg: float, TorRmsAC: float, TorRippleRms: float) -> None: - print("\nTorque mean value: {:.2f} Nm".format(TorAvg)) - # print("\nTorque ripple mean value: {:.2f} Nm".format(TorAvgAC)) - print("\nTorque ripple rms value: {:.2f} Nm".format(TorRmsAC)) - # print("\nTorque ripple relative value: {:.2f} %".format(TorRippleAvg)) - print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) diff --git a/generate_designs.ipynb b/generate_designs.ipynb index 6a2ca34..3cf8201 100644 --- a/generate_designs.ipynb +++ b/generate_designs.ipynb @@ -8,9 +8,8 @@ "source": [ "import os\n", "import numpy as np\n", - "import matplotlib.pyplot as plt\n", "from design import Design\n", - "from geometry import create_barrier\n", + "from geometry import analyze_results, plot_barriers\n", "\n", "project_name = \"SynRM_test\"\n", "design_name = \"Design01\"\n", @@ -19,7 +18,6 @@ "for path in [path_data, path_results]:\n", " os.makedirs(path, exist_ok=True)\n", "file_name_aedt = f'{path_data}/{project_name}.aedt'\n", - "save_design = False\n", "plot_design = True\n", "n_designs = 3\n", "\n", @@ -53,64 +51,33 @@ "metadata": {}, "outputs": [], "source": [ - "w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0\n", + "from generators import FourStupid\n", "\n", - "# Loop over the number of designs\n", + "generator = FourStupid(design)\n", + "\n", + "n_designs = 1\n", "for i in range(0, n_designs):\n", - " # Parameterize the design\n", - " w_mins = w_mins_base + 0.5*np.random.random(4)\n", - " w_mids = w_mins\n", - " y_min0 = 13+4*np.random.random()\n", - " y_min1 = y_min0+w_mins[0]+2+3*np.random.random()\n", - " y_min2 = y_min1+w_mins[1]+2+3*np.random.random()\n", - " y_min3 = y_min2+w_mins[2]+2+2*np.random.random()\n", - " y_mins = np.array([y_min0, y_min1, y_min2, y_min3])\n", - " y_mids = y_mins + np.array([2, 1.5, 1, 0.5])\n", - " thetas = [1, 8, 15, 21]\n", - " w_maxs = w_mins - np.array([0.5, 0.5, 0.5, 0])\n", - " \n", " # Add rotor\n", " design.add_rotor()\n", "\n", " # Add barriers\n", - " for args in zip(y_mins, w_mins, y_mids, w_mids, thetas, w_maxs):\n", - " x_all, y_all = create_barrier(design, *args)\n", - " xy_all = np.vstack((x_all, y_all)).T\n", - " design.add_rotor_barriers(xy_all)\n", + " barriers = generator.random_barriers()\n", + " for barrier in barriers:\n", + " design.add_rotor_barrier(barrier)\n", "\n", " # Compute the torque\n", " Tor = design.compute(NUM_CORES)\n", - " TorAvg, TorRmsAC, TorRippleRms = design.analyze_results(Tor)\n", - "\n", - " # Potentially save the design\n", - " file_design = f'{path_results}/adesign_{i}'\n", - " if save_design:\n", - " design.save_design(file_design)\n", + " TorAvg, _, TorRippleRms = analyze_results(Tor)\n", "\n", " # Delete the rotor\n", " design.delete_rotor()\n", "\n", - " # Potentially plot the desing\n", + " # Potentially save the design\n", " if plot_design:\n", - " plt.figure()\n", - " for args in zip(y_mins, w_mins, y_mids, w_mids, thetas, w_maxs):\n", - " x_all, y_all = create_barrier(design, *args)\n", - " plt.plot(x_all, y_all)\n", - " alphas = np.linspace(0*np.pi,2*np.pi/4,100)\n", - " for r in [design.rotor_r_min, design.rotor_r_max]:\n", - " plt.plot(r*np.cos(alphas), r*np.sin(alphas))\n", - " plt.title(f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %')\n", - " plt.savefig(f'{file_design}.png')\n", - " plt.close()\n", - "\n", - " # Save the points\n", - " np.savez(f'{file_design}.npz',\n", - " y_mins=y_mins,\n", - " w_mins=w_mins,\n", - " y_mids=y_mids,\n", - " w_mids=w_mids,\n", - " thetas=thetas,\n", - " w_maxs=w_maxs)" + " title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %'\n", + " file_name = f'{path_results}/design_{i}'\n", + " plot_barriers(barriers, design, title=title, file_name = f\"{file_name}.png\")\n", + " generator.save_barriers(f\"{file_name}.npz\")" ] }, { diff --git a/generators.py b/generators.py new file mode 100644 index 0000000..c425105 --- /dev/null +++ b/generators.py @@ -0,0 +1,98 @@ +from abc import ABC, abstractmethod +import numpy as np +from scipy.interpolate import CubicSpline +from geometry import rotate + + +class BarrierGenerator(ABC): + @abstractmethod + def random_barriers(self) -> list[np.ndarray]: + pass + + @abstractmethod + def save_barriers(self, file_name: str) -> None: + pass + + +class FourStupid(BarrierGenerator): + def __init__(self, design, n=100, c=0.98, der1=1., der2=1., symmetric=True) -> None: + self.design = design + self.n = n + self.c = c + self.der1 = der1 + self.der2 = der2 + self.w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0 + self.symmetric = symmetric + + def create_barrier( + self, + y_min, + w_min, + y_mid, + w_mid, + theta, + w_max, + ): + + rotor_r_max = self.design.rotor_r_max + + theta1 = (theta+45) / 180*np.pi + x_max1 = self.c*rotor_r_max*np.cos(theta1) + y_max1 = self.c*rotor_r_max*np.sin(theta1) + + x1 = [0, x_max1/2, x_max1] + y1 = [y_min, y_mid, y_max1] + f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,self.der1))) + + theta2 = theta1 + w_max / (2*np.pi*rotor_r_max) * 2*np.pi + x_max2 = self.c*rotor_r_max*np.cos(theta2) + y_max2 = self.c*rotor_r_max*np.sin(theta2) + + s = w_mid / np.sqrt(1 + self.der1**2/4) + x2 = [0, x1[1]-s*self.der1/2, x_max2] + y2 = [y_min+w_min, y1[1]+s, y_max2] + f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,self.der2))) + + x_interp1 = np.linspace(x1[0], x1[-1], self.n) + x_interp2 = np.linspace(x2[0], x2[-1], self.n) + x_all = np.concatenate((x_interp1, x_interp2[::-1])) + y_all = np.concatenate((f1(x_interp1), f2(x_interp2)[::-1])) + if self.symmetric: + x_all = np.concatenate((x_all, -x_all[::-1][1:])) + y_all = np.concatenate((y_all, y_all[::-1][1:])) + x_all, y_all = rotate(x_all, y_all, -45) + return x_all, y_all + + def generate_parameters(self): + # Parameterize the design + self.w_mins = self.w_mins_base + 0.5*np.random.random(4) + y_min0 = 13+4*np.random.random() + y_min1 = y_min0+self.w_mins[0]+2+3*np.random.random() + y_min2 = y_min1+self.w_mins[1]+2+3*np.random.random() + y_min3 = y_min2+self.w_mins[2]+2+2*np.random.random() + self.w_mids = self.w_mins + self.y_mins = np.array([y_min0, y_min1, y_min2, y_min3]) + self.y_mids = self.y_mins + np.array([2, 1.5, 1, 0.5]) + self.thetas = [1, 8, 15, 21] + self.w_maxs = self.w_mins - np.array([0.5, 0.5, 0.5, 0]) + + def random_barriers(self) -> list[np.ndarray]: + # Generate internal parameters + self.generate_parameters() + + # Create the barriers + barriers = [] + for args in zip(self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs): + x_all, y_all = self.create_barrier(*args) + xy_all = np.vstack((x_all, y_all)).T + barriers.append(xy_all) + return barriers + + def save_barriers(self, file_name: str): + np.savez(file_name, + y_mins=self.y_mins, + w_mins=w_mins, + y_mids=y_mids, + w_mids=w_mids, + thetas=thetas, + w_maxs=w_maxs) \ No newline at end of file diff --git a/geometry.py b/geometry.py index dfd4d78..16dc2e9 100644 --- a/geometry.py +++ b/geometry.py @@ -1,5 +1,6 @@ import numpy as np -from scipy.interpolate import CubicSpline +import matplotlib.pyplot as plt + def rotation_matrix(alpha, rad=False): if not rad: @@ -11,44 +12,28 @@ def rotate(x, y, alpha, **kwargs): xy_rot = T @ np.vstack((x, y)) return xy_rot[0], xy_rot[1] -def create_barrier( - design, - y_min, - w_min, - y_mid, - w_mid, - theta, - w_max, - n=100, - c=0.98, - symmetric=True - ): - - theta1 = (theta+45) / 180*np.pi - x_max1 = c*design.rotor_r_max*np.cos(theta1) - y_max1 = c*design.rotor_r_max*np.sin(theta1) - der1 = 1 - - x1 = [0, x_max1/2, x_max1] - y1 = [y_min, y_mid, y_max1] - f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,der1))) +def plot_barriers(barriers, design, title=None, file_name=None): + plt.figure() + alphas = np.linspace(0*np.pi,2*np.pi/4,100) + for r in [design.rotor_r_min, design.rotor_r_max]: + plt.plot(r*np.cos(alphas), r*np.sin(alphas)) + for barrier in barriers: + plt.plot(barrier[:,0], barrier[:,1]) + if title is not None: + plt.title(title) + if file_name is not None: + plt.savefig(file_name) + plt.close() - theta2 = theta1 + w_max / (2*np.pi*design.rotor_r_max) * 2*np.pi - x_max2 = c*design.rotor_r_max*np.cos(theta2) - y_max2 = c*design.rotor_r_max*np.sin(theta2) - der2 = 1 - - s = w_mid / np.sqrt(1 + der1**2/4) - x2 = [0, x1[1]-s*der1/2, x_max2] - y2 = [y_min+w_min, y1[1]+s, y_max2] - f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,der2))) - - x_interp1 = np.linspace(x1[0], x1[-1], n) - x_interp2 = np.linspace(x2[0], x2[-1], n) - x_all = np.concatenate((x_interp1, x_interp2[::-1])) - y_all = np.concatenate((f1(x_interp1), f2(x_interp2)[::-1])) - if symmetric: - x_all = np.concatenate((x_all, -x_all[::-1][1:])) - y_all = np.concatenate((y_all, y_all[::-1][1:])) - x_all, y_all = rotate(x_all, y_all, -45) - return x_all, y_all \ No newline at end of file +def analyze_results(Tor: np.ndarray) -> tuple[float, float, float]: + TorAvg = np.mean(Tor[:-1]) + TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) + TorRippleRms = TorRmsAC/TorAvg*100 + return TorAvg, TorRmsAC, TorRippleRms + +def print_results(TorAvg: float, TorRmsAC: float, TorRippleRms: float) -> None: + print("\nTorque mean value: {:.2f} Nm".format(TorAvg)) + # print("\nTorque ripple mean value: {:.2f} Nm".format(TorAvgAC)) + print("\nTorque ripple rms value: {:.2f} Nm".format(TorRmsAC)) + # print("\nTorque ripple relative value: {:.2f} %".format(TorRippleAvg)) + print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) \ No newline at end of file From c43df9b23570347de2fa119ce94e77c37ca53518 Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 17:12:21 +0100 Subject: [PATCH 06/55] Values fixed --- generators.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/generators.py b/generators.py index c425105..384df6f 100644 --- a/generators.py +++ b/generators.py @@ -91,8 +91,8 @@ def random_barriers(self) -> list[np.ndarray]: def save_barriers(self, file_name: str): np.savez(file_name, y_mins=self.y_mins, - w_mins=w_mins, - y_mids=y_mids, - w_mids=w_mids, - thetas=thetas, - w_maxs=w_maxs) \ No newline at end of file + w_mins=self.w_mins, + y_mids=self.y_mids, + w_mids=self.w_mids, + thetas=self.thetas, + w_maxs=self.w_maxs) \ No newline at end of file From 9212435b4093b68bfb617e139b861777e9c3fdba Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 17:51:58 +0100 Subject: [PATCH 07/55] Added HacklGenerator --- generators.py | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 77 insertions(+), 1 deletion(-) diff --git a/generators.py b/generators.py index 384df6f..20b148e 100644 --- a/generators.py +++ b/generators.py @@ -95,4 +95,80 @@ def save_barriers(self, file_name: str): y_mids=self.y_mids, w_mids=self.w_mids, thetas=self.thetas, - w_maxs=self.w_maxs) \ No newline at end of file + w_maxs=self.w_maxs) + + +class HacklGenerator(BarrierGenerator): + def __init__(self, design): + dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") + airgap = design.mm_to_str("geom_params", "Airgap") + self.R = (dia_stator_gap / 2.0) - airgap - 0.7 + self.phis_inner_min = np.asarray([2., 16., 28.]) + self.phis_inner_max = np.asarray([6., 20., 32.]) + self.phis_outer_min = np.asarray([8., 22., 34.]) + self.phis_outer_max = np.asarray([12., 26., 38.]) + self.n_barriers = len(self.phis_inner_min) + + def generate_parameters(self): + self.lam = 0.25 + self.phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) + self.phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) + + def generate_parameters_representative(self): + self.lam = 0.25 + self.phis_inner = [4.0, 17.3, 30.6] + self.phis_outer = [10.6, 24.0, 37.3] + + def get_arc(self, start_deg, end_deg, num_points=15): + angles = np.linspace(np.radians(start_deg), np.radians(end_deg), num_points) + return [np.array([self.R * np.cos(a), self.R * np.sin(a)]) for a in angles] + + def get_bezier_curve(self, phi_deg, num_points=30): + phi_rad = np.radians(phi_deg) + x0 = self.R * np.cos(phi_rad) + y0 = self.R * np.sin(phi_rad) + r0 = np.array([x0, y0]) + + # Formula: r_{k,1} = [[lam, 1-lam], [0, 1]] * r_{k,0} + x1 = self.lam * x0 + (1 - self.lam) * y0 + y1 = y0 + r1 = np.array([x1, y1]) + + # Formula: r_{k,3} and r_{k,2} (mirrored across the y=x diagonal) + r3 = np.array([y0, x0]) + r2 = np.array([y1, x1]) + + # Formula 33: Polynomial expansion to generate continuous points + z_vals = np.linspace(0, 1, num_points) + curve_pts = [] + for z in z_vals: + pt = (1-z)**3 * r0 + 3*(1-z)**2 * z * r1 + 3*(1-z) * z**2 * r2 + z**3 * r3 + curve_pts.append(pt) + + return curve_pts + + def random_barriers(self) -> list[np.ndarray]: + # Generate internal parameters + self.generate_parameters() + + # Create the barriers + barriers = [] + for phi_inner, phi_outer in zip(self.phis_inner, self.phis_outer): + # Generate the longer part of the barrier + pts_outer = self.get_bezier_curve(phi_outer) + pts_inner = self.get_bezier_curve(phi_inner) + + # Generate the arcs at the end + arc_top = self.get_arc(90 - phi_outer, 90 - phi_inner) + arc_bottom = self.get_arc(phi_inner, phi_outer) + + # Merge them together + closed_loop_pts = pts_outer + arc_top[1:] + pts_inner[::-1][1:] + arc_bottom[1:-1] + barriers.append(np.asarray(closed_loop_pts)) + return barriers + + def save_barriers(self, file_name: str): + np.savez(file_name, + phis_inner=self.phis_inner, + phis_outer=self.phis_outer, + lam=self.lam) From 8407d6e8f035005c44c541ba91cd4e525b26a45d Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 4 Mar 2026 18:20:58 +0100 Subject: [PATCH 08/55] Generation simplified --- generators.py | 46 +++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/generators.py b/generators.py index 20b148e..68d8981 100644 --- a/generators.py +++ b/generators.py @@ -121,31 +121,31 @@ def generate_parameters_representative(self): def get_arc(self, start_deg, end_deg, num_points=15): angles = np.linspace(np.radians(start_deg), np.radians(end_deg), num_points) - return [np.array([self.R * np.cos(a), self.R * np.sin(a)]) for a in angles] + + x = self.R * np.cos(angles) + y = self.R * np.sin(angles) + + return np.column_stack((x, y)) def get_bezier_curve(self, phi_deg, num_points=30): phi_rad = np.radians(phi_deg) - x0 = self.R * np.cos(phi_rad) - y0 = self.R * np.sin(phi_rad) - r0 = np.array([x0, y0]) - - # Formula: r_{k,1} = [[lam, 1-lam], [0, 1]] * r_{k,0} - x1 = self.lam * x0 + (1 - self.lam) * y0 - y1 = y0 - r1 = np.array([x1, y1]) - - # Formula: r_{k,3} and r_{k,2} (mirrored across the y=x diagonal) - r3 = np.array([y0, x0]) - r2 = np.array([y1, x1]) - + x_end = self.R * np.cos(phi_rad) + y_end = self.R * np.sin(phi_rad) + x_bezier = self.lam * x_end + (1 - self.lam) * y_end + + r0 = np.array([x_end, y_end]) + r1 = np.array([x_bezier, y_end]) + r2 = np.array([y_end, x_bezier]) + r3 = np.array([y_end, x_end]) + # Formula 33: Polynomial expansion to generate continuous points - z_vals = np.linspace(0, 1, num_points) - curve_pts = [] - for z in z_vals: - pt = (1-z)**3 * r0 + 3*(1-z)**2 * z * r1 + 3*(1-z) * z**2 * r2 + z**3 * r3 - curve_pts.append(pt) - - return curve_pts + z_vals = np.linspace(0, 1, num_points)[:, None] + return ( + (1 - z_vals)**3 * r0 + + 3 * (1 - z_vals)**2 * z_vals * r1 + + 3 * (1 - z_vals) * z_vals**2 * r2 + + z_vals**3 * r3 + ) def random_barriers(self) -> list[np.ndarray]: # Generate internal parameters @@ -163,8 +163,8 @@ def random_barriers(self) -> list[np.ndarray]: arc_bottom = self.get_arc(phi_inner, phi_outer) # Merge them together - closed_loop_pts = pts_outer + arc_top[1:] + pts_inner[::-1][1:] + arc_bottom[1:-1] - barriers.append(np.asarray(closed_loop_pts)) + barrier = np.concat((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) + barriers.append(barrier) return barriers def save_barriers(self, file_name: str): From 30ea21b1dd2a06e9ec5615b04ca0f07ef0b80ad8 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 09:16:15 +0100 Subject: [PATCH 09/55] Added name property --- generators.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/generators.py b/generators.py index 68d8981..061eab7 100644 --- a/generators.py +++ b/generators.py @@ -13,6 +13,10 @@ def random_barriers(self) -> list[np.ndarray]: def save_barriers(self, file_name: str) -> None: pass + @property + def name(self) -> str: + return self.__class__.__name__ + class FourStupid(BarrierGenerator): def __init__(self, design, n=100, c=0.98, der1=1., der2=1., symmetric=True) -> None: From ac346ddac784ffc570cebde993229992fe03b8b7 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 10:40:19 +0100 Subject: [PATCH 10/55] Added axis equal --- geometry.py | 1 + 1 file changed, 1 insertion(+) diff --git a/geometry.py b/geometry.py index 16dc2e9..328ba07 100644 --- a/geometry.py +++ b/geometry.py @@ -19,6 +19,7 @@ def plot_barriers(barriers, design, title=None, file_name=None): plt.plot(r*np.cos(alphas), r*np.sin(alphas)) for barrier in barriers: plt.plot(barrier[:,0], barrier[:,1]) + plt.axis("equal") if title is not None: plt.title(title) if file_name is not None: From b666f7a51cce8b4a90ba4be1b504f775c3654390 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 10:40:43 +0100 Subject: [PATCH 11/55] Possibility to split the barriers --- generators.py | 117 +++++++++++++++++++++++++++++++++++++------------- 1 file changed, 86 insertions(+), 31 deletions(-) diff --git a/generators.py b/generators.py index 061eab7..6e22e68 100644 --- a/generators.py +++ b/generators.py @@ -4,9 +4,24 @@ from geometry import rotate +def signed_distance(P, a, b, c): + return a * P[:,0] + b * P[:,1] + c + +def check_barrier(barrier: np.ndarray) -> None: + if not np.array_equal(barrier[0], barrier[-1]): + raise ValueError("Barrier must start and end at the same point") + + class BarrierGenerator(ABC): + def __init__(self, offset: float | None = None) -> None: + self.offset = offset + @abstractmethod - def random_barriers(self) -> list[np.ndarray]: + def generate_barriers(self) -> list[np.ndarray]: + pass + + @abstractmethod + def generate_parameters(self) -> None: pass @abstractmethod @@ -17,9 +32,51 @@ def save_barriers(self, file_name: str) -> None: def name(self) -> str: return self.__class__.__name__ + def random_barriers(self) -> list[np.ndarray]: + self.generate_parameters() + return self.generate_barriers() + + def split_barrier(self, barrier: np.ndarray, n_points: int = 20) -> list[np.ndarray]: + assert self.offset is not None + check_barrier(barrier) + + lam = np.linspace(0, 1, n_points)[:, None] + barriers = [] + for a, b, c in zip([1, -1], [-1, 1], [-self.offset, -self.offset]): + f = signed_distance(barrier, a, b, c) + f_is_positive = f > 0 + # TODO: add a better interpolation to the boundary + if f_is_positive[0]: + mask_negative = np.where(~f_is_positive)[0] + idx_min = np.min(mask_negative) + idx_max = np.max(mask_negative) + if len(np.unique(f_is_positive[idx_min:idx_max])) != 1: + raise ValueError("The function crosses the center line too many times.") + barrier_part1 = barrier[:idx_min-1] + barrier_part2 = (1-lam)*barrier[idx_min-1] + lam*barrier[idx_max+1] + barrier_part3 = barrier[idx_max+2:] + barrier_new = np.concatenate((barrier_part1, barrier_part2, barrier_part3)) + else: + mask_positive = np.where(f_is_positive)[0] + idx_min = np.min(mask_positive) + idx_max = np.max(mask_positive) + if len(np.unique(f_is_positive[idx_min:idx_max])) != 1: + raise ValueError("The function crosses the center line too many times.") + barrier_part1 = barrier[idx_min:idx_max] + barrier_part2 = (1-lam)*barrier[idx_max] + lam*barrier[idx_min] + barrier_new = np.concatenate((barrier_part1, barrier_part2)) + check_barrier(barrier_new) + barriers.append(barrier_new) + return barriers + + def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: + if self.offset is None: + return barriers + else: + return [x for barrier in barriers for x in self.split_barrier(barrier)] class FourStupid(BarrierGenerator): - def __init__(self, design, n=100, c=0.98, der1=1., der2=1., symmetric=True) -> None: + def __init__(self, design, n=300, c=0.98, der1=1., der2=1., symmetric=True, **kwargs) -> None: self.design = design self.n = n self.c = c @@ -27,6 +84,7 @@ def __init__(self, design, n=100, c=0.98, der1=1., der2=1., symmetric=True) -> N self.der2 = der2 self.w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0 self.symmetric = symmetric + super().__init__(**kwargs) def create_barrier( self, @@ -67,7 +125,7 @@ def create_barrier( x_all, y_all = rotate(x_all, y_all, -45) return x_all, y_all - def generate_parameters(self): + def generate_parameters(self) -> None: # Parameterize the design self.w_mins = self.w_mins_base + 0.5*np.random.random(4) y_min0 = 13+4*np.random.random() @@ -80,11 +138,7 @@ def generate_parameters(self): self.thetas = [1, 8, 15, 21] self.w_maxs = self.w_mins - np.array([0.5, 0.5, 0.5, 0]) - def random_barriers(self) -> list[np.ndarray]: - # Generate internal parameters - self.generate_parameters() - - # Create the barriers + def generate_barriers(self) -> list[np.ndarray]: barriers = [] for args in zip(self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs): x_all, y_all = self.create_barrier(*args) @@ -92,7 +146,7 @@ def random_barriers(self) -> list[np.ndarray]: barriers.append(xy_all) return barriers - def save_barriers(self, file_name: str): + def save_barriers(self, file_name: str) -> None: np.savez(file_name, y_mins=self.y_mins, w_mins=self.w_mins, @@ -103,7 +157,7 @@ def save_barriers(self, file_name: str): class HacklGenerator(BarrierGenerator): - def __init__(self, design): + def __init__(self, design, **kwargs): dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") airgap = design.mm_to_str("geom_params", "Airgap") self.R = (dia_stator_gap / 2.0) - airgap - 0.7 @@ -112,6 +166,7 @@ def __init__(self, design): self.phis_outer_min = np.asarray([8., 22., 34.]) self.phis_outer_max = np.asarray([12., 26., 38.]) self.n_barriers = len(self.phis_inner_min) + super().__init__(**kwargs) def generate_parameters(self): self.lam = 0.25 @@ -123,6 +178,26 @@ def generate_parameters_representative(self): self.phis_inner = [4.0, 17.3, 30.6] self.phis_outer = [10.6, 24.0, 37.3] + def generate_barriers(self) -> list[np.ndarray]: + # Generate internal parameters + self.generate_parameters() + + # Create the barriers + barriers = [] + for phi_inner, phi_outer in zip(self.phis_inner, self.phis_outer): + # Generate the longer part of the barrier + pts_outer = self.get_bezier_curve(phi_outer) + pts_inner = self.get_bezier_curve(phi_inner) + + # Generate the arcs at the end + arc_top = self.get_arc(90 - phi_outer, 90 - phi_inner) + arc_bottom = self.get_arc(phi_inner, phi_outer) + + # Merge them together + barrier = np.concat((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) + barriers.append(barrier) + return barriers + def get_arc(self, start_deg, end_deg, num_points=15): angles = np.linspace(np.radians(start_deg), np.radians(end_deg), num_points) @@ -131,7 +206,7 @@ def get_arc(self, start_deg, end_deg, num_points=15): return np.column_stack((x, y)) - def get_bezier_curve(self, phi_deg, num_points=30): + def get_bezier_curve(self, phi_deg, num_points=300): phi_rad = np.radians(phi_deg) x_end = self.R * np.cos(phi_rad) y_end = self.R * np.sin(phi_rad) @@ -151,26 +226,6 @@ def get_bezier_curve(self, phi_deg, num_points=30): + z_vals**3 * r3 ) - def random_barriers(self) -> list[np.ndarray]: - # Generate internal parameters - self.generate_parameters() - - # Create the barriers - barriers = [] - for phi_inner, phi_outer in zip(self.phis_inner, self.phis_outer): - # Generate the longer part of the barrier - pts_outer = self.get_bezier_curve(phi_outer) - pts_inner = self.get_bezier_curve(phi_inner) - - # Generate the arcs at the end - arc_top = self.get_arc(90 - phi_outer, 90 - phi_inner) - arc_bottom = self.get_arc(phi_inner, phi_outer) - - # Merge them together - barrier = np.concat((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) - barriers.append(barrier) - return barriers - def save_barriers(self, file_name: str): np.savez(file_name, phis_inner=self.phis_inner, From 65b6f37e9cbe6e5d4bc7bddfb7b83833c41e7326 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 13:06:11 +0000 Subject: [PATCH 12/55] Fixed bounds --- generators.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/generators.py b/generators.py index 6e22e68..bbc911c 100644 --- a/generators.py +++ b/generators.py @@ -161,15 +161,17 @@ def __init__(self, design, **kwargs): dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") airgap = design.mm_to_str("geom_params", "Airgap") self.R = (dia_stator_gap / 2.0) - airgap - 0.7 - self.phis_inner_min = np.asarray([2., 16., 28.]) - self.phis_inner_max = np.asarray([6., 20., 32.]) - self.phis_outer_min = np.asarray([8., 22., 34.]) - self.phis_outer_max = np.asarray([12., 26., 38.]) + self.phis_inner_min = np.asarray([4., 17.3, 30.6]) + self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) + self.phis_outer_min = np.asarray([10.6, 24., 37.3]) + self.phis_outer_max = np.asarray([16.8, 30.2, 34.]) + self.lam_min = 0.25 + self.lam_max = 0.45 self.n_barriers = len(self.phis_inner_min) super().__init__(**kwargs) def generate_parameters(self): - self.lam = 0.25 + self.lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] self.phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) self.phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) @@ -194,7 +196,7 @@ def generate_barriers(self) -> list[np.ndarray]: arc_bottom = self.get_arc(phi_inner, phi_outer) # Merge them together - barrier = np.concat((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) + barrier = np.concatenate((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) barriers.append(barrier) return barriers From 142e237095aec94ae8690a223ca96fe70733ae2b Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 13:19:47 +0000 Subject: [PATCH 13/55] Presribed correct distances --- generators.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/generators.py b/generators.py index bbc911c..1924310 100644 --- a/generators.py +++ b/generators.py @@ -11,6 +11,11 @@ def check_barrier(barrier: np.ndarray) -> None: if not np.array_equal(barrier[0], barrier[-1]): raise ValueError("Barrier must start and end at the same point") +def compute_r_max(design, r_stator_end: float) -> float: + dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") + airgap = design.mm_to_str("geom_params", "Airgap") + return (dia_stator_gap / 2.0) - airgap - r_stator_end + class BarrierGenerator(ABC): def __init__(self, offset: float | None = None) -> None: @@ -76,10 +81,10 @@ def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: return [x for barrier in barriers for x in self.split_barrier(barrier)] class FourStupid(BarrierGenerator): - def __init__(self, design, n=300, c=0.98, der1=1., der2=1., symmetric=True, **kwargs) -> None: + def __init__(self, design, r_stator_end, n=300, der1=1., der2=1., symmetric=True, **kwargs) -> None: + self.R = compute_r_max(design, r_stator_end) self.design = design self.n = n - self.c = c self.der1 = der1 self.der2 = der2 self.w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0 @@ -96,19 +101,17 @@ def create_barrier( w_max, ): - rotor_r_max = self.design.rotor_r_max - theta1 = (theta+45) / 180*np.pi - x_max1 = self.c*rotor_r_max*np.cos(theta1) - y_max1 = self.c*rotor_r_max*np.sin(theta1) + x_max1 = self.R*np.cos(theta1) + y_max1 = self.R*np.sin(theta1) x1 = [0, x_max1/2, x_max1] y1 = [y_min, y_mid, y_max1] f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,self.der1))) - theta2 = theta1 + w_max / (2*np.pi*rotor_r_max) * 2*np.pi - x_max2 = self.c*rotor_r_max*np.cos(theta2) - y_max2 = self.c*rotor_r_max*np.sin(theta2) + theta2 = theta1 + w_max / self.design.rotor_r_max + x_max2 = self.R*np.cos(theta2) + y_max2 = self.R*np.sin(theta2) s = w_mid / np.sqrt(1 + self.der1**2/4) x2 = [0, x1[1]-s*self.der1/2, x_max2] @@ -157,10 +160,8 @@ def save_barriers(self, file_name: str) -> None: class HacklGenerator(BarrierGenerator): - def __init__(self, design, **kwargs): - dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") - airgap = design.mm_to_str("geom_params", "Airgap") - self.R = (dia_stator_gap / 2.0) - airgap - 0.7 + def __init__(self, design, r_stator_end, **kwargs): + self.R = compute_r_max(design, r_stator_end) self.phis_inner_min = np.asarray([4., 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) self.phis_outer_min = np.asarray([10.6, 24., 37.3]) From 7f4e333ebf947824e757b770369032656ba21f8d Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 13:30:57 +0000 Subject: [PATCH 14/55] Clean up --- example.py | 45 --------------------------------------------- 1 file changed, 45 deletions(-) delete mode 100644 example.py diff --git a/example.py b/example.py deleted file mode 100644 index 04db18f..0000000 --- a/example.py +++ /dev/null @@ -1,45 +0,0 @@ -import os -from ansys.aedt.core import Maxwell2d -from design import Design - -project_name = "SynRM_test" -design_name = "Design01" -path = os.path.join(os.getcwd(), 'data') -os.makedirs(path, exist_ok=True) -file_name_aedt = f'{path}/{project_name}.aedt' - -# Define constants -AEDT_VERSION = "2024.2" -NUM_CORES = 4 -NG_MODE = True #non-graphical mode -CLS_EXIT = True #close on exit - -barrier_points = [ - [14, 14], - [34, 2], - [36, 2], - [20, 20], - [2, 36], - [2, 34], - [14, 14], -] - -if not os.path.exists(file_name_aedt): - m2d = Maxwell2d(project=project_name, design=design_name, version=AEDT_VERSION, - non_graphical=NG_MODE, new_desktop=False, - close_on_exit=CLS_EXIT, student_version=False, solution_type="TransientXY") - design = Design(m2d) - design.create_stator() - design.save_project(file_name_aedt) -else: - design = Design() - design.load_stator(file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, - new_desktop=False, close_on_exit=CLS_EXIT, student_version=False) - -design.add_rotor() -design.add_rotor_holes(barrier_points) -Tor = design.compute(NUM_CORES) -results = design.analyze_results(Tor) -design.print_results(*results) -design.delete_rotor() -design.close_project() From 0dfd6a4aa4d72f8cbb7f2839eca508041f0a4db3 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 13:35:25 +0000 Subject: [PATCH 15/55] Added typing --- generators.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/generators.py b/generators.py index 1924310..022cca8 100644 --- a/generators.py +++ b/generators.py @@ -1,24 +1,25 @@ from abc import ABC, abstractmethod import numpy as np from scipy.interpolate import CubicSpline +from design import Design from geometry import rotate -def signed_distance(P, a, b, c): +def signed_distance(P: np.ndarray, a: float, b: float, c: float) -> np.ndarray: return a * P[:,0] + b * P[:,1] + c def check_barrier(barrier: np.ndarray) -> None: if not np.array_equal(barrier[0], barrier[-1]): raise ValueError("Barrier must start and end at the same point") -def compute_r_max(design, r_stator_end: float) -> float: +def compute_r_max(design: Design, r_stator_end: float) -> float: dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") airgap = design.mm_to_str("geom_params", "Airgap") return (dia_stator_gap / 2.0) - airgap - r_stator_end class BarrierGenerator(ABC): - def __init__(self, offset: float | None = None) -> None: + def __init__(self, offset: float) -> None: self.offset = offset @abstractmethod @@ -47,7 +48,7 @@ def split_barrier(self, barrier: np.ndarray, n_points: int = 20) -> list[np.ndar lam = np.linspace(0, 1, n_points)[:, None] barriers = [] - for a, b, c in zip([1, -1], [-1, 1], [-self.offset, -self.offset]): + for a, b, c in zip([1., -1.], [-1., 1.], [-self.offset, -self.offset]): f = signed_distance(barrier, a, b, c) f_is_positive = f > 0 # TODO: add a better interpolation to the boundary @@ -75,10 +76,7 @@ def split_barrier(self, barrier: np.ndarray, n_points: int = 20) -> list[np.ndar return barriers def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: - if self.offset is None: - return barriers - else: - return [x for barrier in barriers for x in self.split_barrier(barrier)] + return [x for barrier in barriers for x in self.split_barrier(barrier)] class FourStupid(BarrierGenerator): def __init__(self, design, r_stator_end, n=300, der1=1., der2=1., symmetric=True, **kwargs) -> None: From dbd49c3766f7bb6642b5779181c5cff04d32a8f9 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 13:40:56 +0000 Subject: [PATCH 16/55] Refactor --- generators.py | 35 ++++++++++++++++------------------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/generators.py b/generators.py index 022cca8..655211d 100644 --- a/generators.py +++ b/generators.py @@ -5,6 +5,14 @@ from geometry import rotate +def get_arc(R: float, start_deg: float, end_deg: float, num_points: int = 15) -> np.ndarray: + angles = np.linspace(np.radians(start_deg), np.radians(end_deg), num_points) + + x = R * np.cos(angles) + y = R * np.sin(angles) + + return np.column_stack((x, y)) + def signed_distance(P: np.ndarray, a: float, b: float, c: float) -> np.ndarray: return a * P[:,0] + b * P[:,1] + c @@ -89,7 +97,7 @@ def __init__(self, design, r_stator_end, n=300, der1=1., der2=1., symmetric=True self.symmetric = symmetric super().__init__(**kwargs) - def create_barrier( + def _create_barrier( self, y_min, w_min, @@ -116,6 +124,7 @@ def create_barrier( y2 = [y_min+w_min, y1[1]+s, y_max2] f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,self.der2))) + # TODO: use get_arc x_interp1 = np.linspace(x1[0], x1[-1], self.n) x_interp2 = np.linspace(x2[0], x2[-1], self.n) x_all = np.concatenate((x_interp1, x_interp2[::-1])) @@ -142,7 +151,7 @@ def generate_parameters(self) -> None: def generate_barriers(self) -> list[np.ndarray]: barriers = [] for args in zip(self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs): - x_all, y_all = self.create_barrier(*args) + x_all, y_all = self._create_barrier(*args) xy_all = np.vstack((x_all, y_all)).T barriers.append(xy_all) return barriers @@ -180,34 +189,22 @@ def generate_parameters_representative(self): self.phis_outer = [10.6, 24.0, 37.3] def generate_barriers(self) -> list[np.ndarray]: - # Generate internal parameters - self.generate_parameters() - - # Create the barriers barriers = [] for phi_inner, phi_outer in zip(self.phis_inner, self.phis_outer): # Generate the longer part of the barrier - pts_outer = self.get_bezier_curve(phi_outer) - pts_inner = self.get_bezier_curve(phi_inner) + pts_outer = self._get_bezier_curve(phi_outer) + pts_inner = self._get_bezier_curve(phi_inner) # Generate the arcs at the end - arc_top = self.get_arc(90 - phi_outer, 90 - phi_inner) - arc_bottom = self.get_arc(phi_inner, phi_outer) + arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner) + arc_bottom = get_arc(self.R, phi_inner, phi_outer) # Merge them together barrier = np.concatenate((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) barriers.append(barrier) return barriers - def get_arc(self, start_deg, end_deg, num_points=15): - angles = np.linspace(np.radians(start_deg), np.radians(end_deg), num_points) - - x = self.R * np.cos(angles) - y = self.R * np.sin(angles) - - return np.column_stack((x, y)) - - def get_bezier_curve(self, phi_deg, num_points=300): + def _get_bezier_curve(self, phi_deg, num_points=300): phi_rad = np.radians(phi_deg) x_end = self.R * np.cos(phi_rad) y_end = self.R * np.sin(phi_rad) From 8eb72301b23faff5f363ba20c8158a05db9b50e0 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 13:44:26 +0000 Subject: [PATCH 17/55] Notebook updated --- generate_designs.ipynb | 49 ++++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 23 deletions(-) diff --git a/generate_designs.ipynb b/generate_designs.ipynb index 3cf8201..7db8ebd 100644 --- a/generate_designs.ipynb +++ b/generate_designs.ipynb @@ -9,6 +9,7 @@ "import os\n", "import numpy as np\n", "from design import Design\n", + "from generators import FourStupid, HacklGenerator\n", "from geometry import analyze_results, plot_barriers\n", "\n", "project_name = \"SynRM_test\"\n", @@ -19,10 +20,10 @@ " os.makedirs(path, exist_ok=True)\n", "file_name_aedt = f'{path_data}/{project_name}.aedt'\n", "plot_design = True\n", - "n_designs = 3\n", + "n_designs = 10\n", "\n", "# Define constants\n", - "AEDT_VERSION = \"2024.2\"\n", + "AEDT_VERSION = \"2024.1\"\n", "NUM_CORES = 4\n", "NG_MODE = True #non-graphical mode\n", "CLS_EXIT = True #close on exit\n", @@ -51,33 +52,35 @@ "metadata": {}, "outputs": [], "source": [ - "from generators import FourStupid\n", + "r_stator_end = 0.7\n", + "offset = 0.7 / 2\n", + "generator_stupid = FourStupid(design, r_stator_end, offset=offset)\n", + "generator_hackl = HacklGenerator(design, r_stator_end, offset=offset)\n", + "generators = [generator_stupid, generator_hackl]\n", "\n", - "generator = FourStupid(design)\n", - "\n", - "n_designs = 1\n", "for i in range(0, n_designs):\n", - " # Add rotor\n", - " design.add_rotor()\n", + " for generator in generators:\n", + " barriers1 = generator.random_barriers()\n", + " barriers2 = generator.split_barriers(barriers1)\n", + " for j, barriers in enumerate([barriers1, barriers2]):\n", + " design.add_rotor()\n", "\n", - " # Add barriers\n", - " barriers = generator.random_barriers()\n", - " for barrier in barriers:\n", - " design.add_rotor_barrier(barrier)\n", + " for barrier in barriers:\n", + " design.add_rotor_barrier(barrier)\n", "\n", - " # Compute the torque\n", - " Tor = design.compute(NUM_CORES)\n", - " TorAvg, _, TorRippleRms = analyze_results(Tor)\n", + " # Compute the torque\n", + " Tor = design.compute(NUM_CORES)\n", + " TorAvg, _, TorRippleRms = analyze_results(Tor)\n", "\n", - " # Delete the rotor\n", - " design.delete_rotor()\n", + " # Delete the rotor\n", + " design.delete_rotor()\n", "\n", - " # Potentially save the design\n", - " if plot_design:\n", - " title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %'\n", - " file_name = f'{path_results}/design_{i}'\n", - " plot_barriers(barriers, design, title=title, file_name = f\"{file_name}.png\")\n", - " generator.save_barriers(f\"{file_name}.npz\")" + " # Potentially save the design\n", + " if plot_design:\n", + " title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %'\n", + " file_name = f'{path_results}/design_{generator.name}_{i}_{j}'\n", + " plot_barriers(barriers, design, title=title, file_name = f\"{file_name}.png\")\n", + " generator.save_barriers(f\"{file_name}.npz\")" ] }, { From 6d03db8db05c07ef94d9fe0a201fec2437d96760 Mon Sep 17 00:00:00 2001 From: sadda Date: Thu, 5 Mar 2026 14:32:53 +0000 Subject: [PATCH 18/55] Fix bounds --- generators.py | 44 +++++++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/generators.py b/generators.py index 655211d..db44957 100644 --- a/generators.py +++ b/generators.py @@ -59,28 +59,26 @@ def split_barrier(self, barrier: np.ndarray, n_points: int = 20) -> list[np.ndar for a, b, c in zip([1., -1.], [-1., 1.], [-self.offset, -self.offset]): f = signed_distance(barrier, a, b, c) f_is_positive = f > 0 - # TODO: add a better interpolation to the boundary - if f_is_positive[0]: - mask_negative = np.where(~f_is_positive)[0] - idx_min = np.min(mask_negative) - idx_max = np.max(mask_negative) - if len(np.unique(f_is_positive[idx_min:idx_max])) != 1: - raise ValueError("The function crosses the center line too many times.") - barrier_part1 = barrier[:idx_min-1] - barrier_part2 = (1-lam)*barrier[idx_min-1] + lam*barrier[idx_max+1] - barrier_part3 = barrier[idx_max+2:] - barrier_new = np.concatenate((barrier_part1, barrier_part2, barrier_part3)) - else: - mask_positive = np.where(f_is_positive)[0] - idx_min = np.min(mask_positive) - idx_max = np.max(mask_positive) - if len(np.unique(f_is_positive[idx_min:idx_max])) != 1: - raise ValueError("The function crosses the center line too many times.") - barrier_part1 = barrier[idx_min:idx_max] - barrier_part2 = (1-lam)*barrier[idx_max] + lam*barrier[idx_min] - barrier_new = np.concatenate((barrier_part1, barrier_part2)) - check_barrier(barrier_new) - barriers.append(barrier_new) + mask_positive = np.where(f_is_positive)[0] + + barrier_new = [barrier[mask_positive[0]]] + for i in range(len(mask_positive) - 1): + i1 = mask_positive[i] + i2 = mask_positive[i+1] + + if i2 - i1 == 1: + barrier_new.append(barrier[i2]) + else: + interp = (1 - lam) * barrier[i1] + lam * barrier[i2] + barrier_new.extend(interp[1:]) + if not np.array_equal(barrier_new[-1], barrier_new[0]): + interp = (1 - lam) * barrier_new[-1] + lam * barrier_new[0] + barrier_new.extend(interp[1:]) + + if len(barrier_new) > 1: + barrier_new = np.array(barrier_new) + check_barrier(barrier_new) + barriers.append(barrier_new) return barriers def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: @@ -172,7 +170,7 @@ def __init__(self, design, r_stator_end, **kwargs): self.phis_inner_min = np.asarray([4., 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) self.phis_outer_min = np.asarray([10.6, 24., 37.3]) - self.phis_outer_max = np.asarray([16.8, 30.2, 34.]) + self.phis_outer_max = np.asarray([16.8, 30.2, 44.]) self.lam_min = 0.25 self.lam_max = 0.45 self.n_barriers = len(self.phis_inner_min) From a5fad0593c26d046cebc7bfeafcac03f61fa0b4d Mon Sep 17 00:00:00 2001 From: sadda Date: Sat, 7 Mar 2026 08:59:19 +0000 Subject: [PATCH 19/55] Gitignore updated --- .gitignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index bf911bb..218a65e 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ data/ -results/ +results*/ +nohup.out # Byte-compiled / optimized / DLL files __pycache__/ From 2781fd43eef3eb003ff91bd8e7d068ce58f4d9c3 Mon Sep 17 00:00:00 2001 From: sadda Date: Sat, 7 Mar 2026 09:02:07 +0000 Subject: [PATCH 20/55] Added AbstractHacklGenerator --- generators.py | 52 ++++++++++++++++++++++++++++++++++----------------- 1 file changed, 35 insertions(+), 17 deletions(-) diff --git a/generators.py b/generators.py index db44957..60c4ec7 100644 --- a/generators.py +++ b/generators.py @@ -164,34 +164,34 @@ def save_barriers(self, file_name: str) -> None: w_maxs=self.w_maxs) -class HacklGenerator(BarrierGenerator): +class AbstractHacklGenerator(BarrierGenerator): def __init__(self, design, r_stator_end, **kwargs): self.R = compute_r_max(design, r_stator_end) self.phis_inner_min = np.asarray([4., 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) self.phis_outer_min = np.asarray([10.6, 24., 37.3]) self.phis_outer_max = np.asarray([16.8, 30.2, 44.]) - self.lam_min = 0.25 - self.lam_max = 0.45 self.n_barriers = len(self.phis_inner_min) super().__init__(**kwargs) - def generate_parameters(self): - self.lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] + @abstractmethod + def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: + pass + + @abstractmethod + def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: + pass + + def generate_parameters(self) -> None: self.phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) self.phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) - def generate_parameters_representative(self): - self.lam = 0.25 - self.phis_inner = [4.0, 17.3, 30.6] - self.phis_outer = [10.6, 24.0, 37.3] - def generate_barriers(self) -> list[np.ndarray]: barriers = [] for phi_inner, phi_outer in zip(self.phis_inner, self.phis_outer): # Generate the longer part of the barrier - pts_outer = self._get_bezier_curve(phi_outer) - pts_inner = self._get_bezier_curve(phi_inner) + pts_outer = self._get_bezier_curve(phi_outer, False) + pts_inner = self._get_bezier_curve(phi_inner, True) # Generate the arcs at the end arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner) @@ -202,15 +202,16 @@ def generate_barriers(self) -> list[np.ndarray]: barriers.append(barrier) return barriers - def _get_bezier_curve(self, phi_deg, num_points=300): + def _get_bezier_curve(self, phi_deg, is_inner, num_points=300): phi_rad = np.radians(phi_deg) x_end = self.R * np.cos(phi_rad) y_end = self.R * np.sin(phi_rad) - x_bezier = self.lam * x_end + (1 - self.lam) * y_end + f_bezier = self._inner_bezier_point if is_inner else self._outer_bezier_point + x_bezier, y_bezier = f_bezier(x_end, y_end) r0 = np.array([x_end, y_end]) - r1 = np.array([x_bezier, y_end]) - r2 = np.array([y_end, x_bezier]) + r1 = np.array([x_bezier, y_bezier]) + r2 = np.array([y_bezier, x_bezier]) r3 = np.array([y_end, x_end]) # Formula 33: Polynomial expansion to generate continuous points @@ -222,8 +223,25 @@ def _get_bezier_curve(self, phi_deg, num_points=300): + z_vals**3 * r3 ) + +class HacklGenerator_OneLambda(AbstractHacklGenerator): + def __init__(self, design, r_stator_end, **kwargs): + self.lam_min = 0.25 + self.lam_max = 0.45 + super().__init__(design, r_stator_end, **kwargs) + + def generate_parameters(self): + self.lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] + super().generate_parameters() + + def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: + return self.lam * x + (1 - self.lam) * y, y + + def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: + return self.lam * x + (1 - self.lam) * y, y + def save_barriers(self, file_name: str): np.savez(file_name, phis_inner=self.phis_inner, phis_outer=self.phis_outer, - lam=self.lam) + lam=self.lam) \ No newline at end of file From 5b64536f87aa1b9fcb69aeafac97d380a1782595 Mon Sep 17 00:00:00 2001 From: sadda Date: Sat, 7 Mar 2026 09:12:36 +0000 Subject: [PATCH 21/55] Handling of n_points --- generators.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/generators.py b/generators.py index 60c4ec7..9c3402d 100644 --- a/generators.py +++ b/generators.py @@ -5,8 +5,8 @@ from geometry import rotate -def get_arc(R: float, start_deg: float, end_deg: float, num_points: int = 15) -> np.ndarray: - angles = np.linspace(np.radians(start_deg), np.radians(end_deg), num_points) +def get_arc(R: float, start_deg: float, end_deg: float, n_points: int) -> np.ndarray: + angles = np.linspace(np.radians(start_deg), np.radians(end_deg), n_points) x = R * np.cos(angles) y = R * np.sin(angles) @@ -27,8 +27,10 @@ def compute_r_max(design: Design, r_stator_end: float) -> float: class BarrierGenerator(ABC): - def __init__(self, offset: float) -> None: + def __init__(self, offset: float, n_curve: int = 500, n_flat: int = 20) -> None: self.offset = offset + self.n_curve = n_curve + self.n_flat = n_flat @abstractmethod def generate_barriers(self) -> list[np.ndarray]: @@ -50,11 +52,11 @@ def random_barriers(self) -> list[np.ndarray]: self.generate_parameters() return self.generate_barriers() - def split_barrier(self, barrier: np.ndarray, n_points: int = 20) -> list[np.ndarray]: + def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: assert self.offset is not None check_barrier(barrier) - lam = np.linspace(0, 1, n_points)[:, None] + lam = np.linspace(0, 1, self.n_flat)[:, None] barriers = [] for a, b, c in zip([1., -1.], [-1., 1.], [-self.offset, -self.offset]): f = signed_distance(barrier, a, b, c) @@ -85,10 +87,9 @@ def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: return [x for barrier in barriers for x in self.split_barrier(barrier)] class FourStupid(BarrierGenerator): - def __init__(self, design, r_stator_end, n=300, der1=1., der2=1., symmetric=True, **kwargs) -> None: + def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwargs) -> None: self.R = compute_r_max(design, r_stator_end) self.design = design - self.n = n self.der1 = der1 self.der2 = der2 self.w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0 @@ -122,9 +123,9 @@ def _create_barrier( y2 = [y_min+w_min, y1[1]+s, y_max2] f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,self.der2))) - # TODO: use get_arc - x_interp1 = np.linspace(x1[0], x1[-1], self.n) - x_interp2 = np.linspace(x2[0], x2[-1], self.n) + x_interp1 = np.linspace(x1[0], x1[-1], self.n_curve) + x_interp2 = np.linspace(x2[0], x2[-1], self.n_curve) + # TODO: use more points for connecting? x_all = np.concatenate((x_interp1, x_interp2[::-1])) y_all = np.concatenate((f1(x_interp1), f2(x_interp2)[::-1])) if self.symmetric: @@ -194,15 +195,15 @@ def generate_barriers(self) -> list[np.ndarray]: pts_inner = self._get_bezier_curve(phi_inner, True) # Generate the arcs at the end - arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner) - arc_bottom = get_arc(self.R, phi_inner, phi_outer) + arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner, self.n_flat) + arc_bottom = get_arc(self.R, phi_inner, phi_outer, self.n_flat) # Merge them together barrier = np.concatenate((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) barriers.append(barrier) return barriers - def _get_bezier_curve(self, phi_deg, is_inner, num_points=300): + def _get_bezier_curve(self, phi_deg, is_inner): phi_rad = np.radians(phi_deg) x_end = self.R * np.cos(phi_rad) y_end = self.R * np.sin(phi_rad) @@ -215,7 +216,7 @@ def _get_bezier_curve(self, phi_deg, is_inner, num_points=300): r3 = np.array([y_end, x_end]) # Formula 33: Polynomial expansion to generate continuous points - z_vals = np.linspace(0, 1, num_points)[:, None] + z_vals = np.linspace(0, 1, self.n_curve)[:, None] return ( (1 - z_vals)**3 * r0 + 3 * (1 - z_vals)**2 * z_vals * r1 From 8b5fb2704ba19ea38eb7b1419fc56f34942420b6 Mon Sep 17 00:00:00 2001 From: sadda Date: Sat, 7 Mar 2026 11:37:56 +0000 Subject: [PATCH 22/55] Added new classes --- generators.py | 57 +++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/generators.py b/generators.py index 9c3402d..aebd0f9 100644 --- a/generators.py +++ b/generators.py @@ -239,10 +239,63 @@ def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: - return self.lam * x + (1 - self.lam) * y, y + return self._inner_bezier_point(x, y) + + def save_barriers(self, file_name: str): + np.savez(file_name, + phis_inner=self.phis_inner, + phis_outer=self.phis_outer, + lam=self.lam) + + + +class HacklGenerator_OneLambdaTheta(HacklGenerator_OneLambda): + def __init__(self, design, r_stator_end, **kwargs): + self.theta_min = -20 + self.theta_max = 20 + super().__init__(design, r_stator_end, **kwargs) + + def generate_parameters(self): + self.theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] + super().generate_parameters() + + def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: + x, y = super()._inner_bezier_point(x, y) + xx, yy = rotate(x, y, self.theta) + return xx[0], yy[0] + + def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: + return self._inner_bezier_point(x, y) + + def save_barriers(self, file_name: str): + np.savez(file_name, + phis_inner=self.phis_inner, + phis_outer=self.phis_outer, + lam=self.lam) + + +class HacklGenerator_TwoLambdas(AbstractHacklGenerator): + def __init__(self, design, r_stator_end, **kwargs): + self.lam_inner_min = 0.25 + self.lam_inner_max = 0.45 + self.lam_outer_min = 0.25 + self.lam_outer_max = 0.45 + super().__init__(design, r_stator_end, **kwargs) + + def generate_parameters(self): + self.lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] + self.lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] + super().generate_parameters() + + def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: + return self.lam_inner * x + (1 - self.lam_inner) * y, y + + def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: + return self.lam_outer * x + (1 - self.lam_outer) * y, y def save_barriers(self, file_name: str): np.savez(file_name, phis_inner=self.phis_inner, phis_outer=self.phis_outer, - lam=self.lam) \ No newline at end of file + lam_inner=self.lam_inner, + lam_outer=self.lam_outer) \ No newline at end of file From 29ca43931fe338ebd27757b787b09a4ad42cf417 Mon Sep 17 00:00:00 2001 From: sadda Date: Sat, 7 Mar 2026 12:09:15 +0000 Subject: [PATCH 23/55] Added more freedom --- generators.py | 46 ++++++++++++++++++---------------------------- 1 file changed, 18 insertions(+), 28 deletions(-) diff --git a/generators.py b/generators.py index aebd0f9..cf6e4e5 100644 --- a/generators.py +++ b/generators.py @@ -176,11 +176,7 @@ def __init__(self, design, r_stator_end, **kwargs): super().__init__(**kwargs) @abstractmethod - def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: - pass - - @abstractmethod - def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: pass def generate_parameters(self) -> None: @@ -189,10 +185,10 @@ def generate_parameters(self) -> None: def generate_barriers(self) -> list[np.ndarray]: barriers = [] - for phi_inner, phi_outer in zip(self.phis_inner, self.phis_outer): + for i, (phi_inner, phi_outer) in enumerate(zip(self.phis_inner, self.phis_outer)): # Generate the longer part of the barrier - pts_outer = self._get_bezier_curve(phi_outer, False) - pts_inner = self._get_bezier_curve(phi_inner, True) + pts_outer = self._get_bezier_curve(phi_outer, False, i) + pts_inner = self._get_bezier_curve(phi_inner, True, i) # Generate the arcs at the end arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner, self.n_flat) @@ -203,12 +199,11 @@ def generate_barriers(self) -> list[np.ndarray]: barriers.append(barrier) return barriers - def _get_bezier_curve(self, phi_deg, is_inner): + def _get_bezier_curve(self, phi_deg, is_inner, i): phi_rad = np.radians(phi_deg) x_end = self.R * np.cos(phi_rad) y_end = self.R * np.sin(phi_rad) - f_bezier = self._inner_bezier_point if is_inner else self._outer_bezier_point - x_bezier, y_bezier = f_bezier(x_end, y_end) + x_bezier, y_bezier = self._bezier_point(x_end, y_end, is_inner, i) r0 = np.array([x_end, y_end]) r1 = np.array([x_bezier, y_bezier]) @@ -235,12 +230,9 @@ def generate_parameters(self): self.lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] super().generate_parameters() - def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y - def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: - return self._inner_bezier_point(x, y) - def save_barriers(self, file_name: str): np.savez(file_name, phis_inner=self.phis_inner, @@ -259,19 +251,19 @@ def generate_parameters(self): self.theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] super().generate_parameters() - def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: - x, y = super()._inner_bezier_point(x, y) - xx, yy = rotate(x, y, self.theta) - return xx[0], yy[0] - - def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: - return self._inner_bezier_point(x, y) + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: + x, y = super()._bezier_point(x, y, is_inner, order) + if order < 2: + xx, yy = rotate(x, y, self.theta) + return xx[0], yy[0] + return x, y def save_barriers(self, file_name: str): np.savez(file_name, phis_inner=self.phis_inner, phis_outer=self.phis_outer, - lam=self.lam) + lam=self.lam, + theta=self.theta) class HacklGenerator_TwoLambdas(AbstractHacklGenerator): @@ -287,11 +279,9 @@ def generate_parameters(self): self.lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] super().generate_parameters() - def _inner_bezier_point(self, x: float, y: float) -> tuple[float, float]: - return self.lam_inner * x + (1 - self.lam_inner) * y, y - - def _outer_bezier_point(self, x: float, y: float) -> tuple[float, float]: - return self.lam_outer * x + (1 - self.lam_outer) * y, y + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: + lam = self.lam_inner if is_inner else self.lam_outer + return lam * x + (1 - lam) * y, y def save_barriers(self, file_name: str): np.savez(file_name, From ac14e8031b4cd4be770de0a13ef14f6c1fc9aca3 Mon Sep 17 00:00:00 2001 From: sadda Date: Sat, 7 Mar 2026 18:19:35 +0000 Subject: [PATCH 24/55] Handling errors --- design.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/design.py b/design.py index 5913dc2..819e974 100644 --- a/design.py +++ b/design.py @@ -508,7 +508,10 @@ def compute(self, NUM_CORES: int = 1): solutions = m2d.post.get_solution_data( expressions="Moving1.Torque", primary_sweep_variable="Time" ) - return solutions.data_magnitude() + try: + return solutions.data_magnitude() + except AttributeError: + return None def delete_rotor(self) -> None: assert isinstance(self.m2d.modeler, Modeler2D) From 6cd3bc58a2e98a9923a81714c4dd417d64a50fb8 Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 10:12:19 +0000 Subject: [PATCH 25/55] Added script --- run.py | 88 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 run.py diff --git a/run.py b/run.py new file mode 100644 index 0000000..c89b15b --- /dev/null +++ b/run.py @@ -0,0 +1,88 @@ +# export ANSYSEM_ROOT241=/data/AnsysEM/v241/Linux64 + +import os +import numpy as np +import pandas as pd +from design import Design +from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta +from geometry import analyze_results, plot_barriers + +project_name = "SynRM_test" +design_name = "Design01" +path_data = os.path.join(os.getcwd(), 'data') +path_results = 'results' +for path in [path_data, path_results]: + os.makedirs(path, exist_ok=True) +file_name_aedt = f'{path_data}/{project_name}.aedt' +plot_design = True +n_designs = 100 + +# Define constants +AEDT_VERSION = "2024.1" +NUM_CORES = 4 +NG_MODE = True #non-graphical mode +CLS_EXIT = True #close on exit + +if not os.path.exists(file_name_aedt): + design = Design.create( + project_name, design_name, file_name_aedt, + version=AEDT_VERSION, + non_graphical=NG_MODE, + new_desktop=False, + close_on_exit=CLS_EXIT, + ) +else: + design = Design.load( + file_name_aedt, + version=AEDT_VERSION, + non_graphical=NG_MODE, + new_desktop=False, + close_on_exit=CLS_EXIT, + ) + +r_stator_end = 0.7 +offset = 0.7 / 2 +generator_stupid = FourStupid(design, r_stator_end, offset=offset) +generator_hackl1 = HacklGenerator_OneLambda(design, r_stator_end, offset=offset) +generator_hackl2 = HacklGenerator_TwoLambdas(design, r_stator_end, offset=offset) +generator_hackl3 = HacklGenerator_OneLambdaTheta(design, r_stator_end, offset=offset) +generators = [generator_stupid, generator_hackl1, generator_hackl2, generator_hackl3] + +metadata = pd.DataFrame() +for i in range(0, n_designs): + for generator in generators: + barriers = generator.random_barriers() + barriers = generator.split_barriers(barriers) + + design.add_rotor() + + for barrier in barriers: + design.add_rotor_barrier(barrier) + + # Compute the torque + Tor = design.compute(NUM_CORES) + if Tor is None: + TorAvg, TorRippleRms = np.nan, np.nan + else: + TorAvg, _, TorRippleRms = analyze_results(Tor) + + # Delete the rotor + design.delete_rotor() + + # Potentially save the design + if plot_design: + title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %' + file_name = f'{path_results}/design_{generator.name}_{i}' + plot_barriers(barriers, design, title=title, file_name = f"{file_name}.png") + generator.save_barriers(f"{file_name}.npz") + + metadata_new = { + 'method': generator.name, + 'design': i, + 'T': TorAvg, + 'ripple': TorRippleRms, + } + metadata = pd.concat((metadata, pd.DataFrame([metadata_new])), ignore_index=True) + metadata.to_csv(f"{path_results}/metadata.csv", index=False) + +design.close_project() \ No newline at end of file From 760a57028e9006a174f37945d7f30a7dff0c9aa8 Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 10:12:27 +0000 Subject: [PATCH 26/55] Bounds exported --- generators.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/generators.py b/generators.py index cf6e4e5..f09aac7 100644 --- a/generators.py +++ b/generators.py @@ -226,6 +226,12 @@ def __init__(self, design, r_stator_end, **kwargs): self.lam_max = 0.45 super().__init__(design, r_stator_end, **kwargs) + @property + def bounds(self) -> tuple[np.ndarray, np.ndarray]: + lb = np.concatenate((self.phis_inner_min, [self.lam_min])) + ub = np.concatenate((self.phis_inner_max, [self.lam_max])) + return lb, ub + def generate_parameters(self): self.lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] super().generate_parameters() From 6bae2b903a77246821a9971950abf26ad6b25886 Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 11:50:29 +0000 Subject: [PATCH 27/55] Refactor --- generators.py | 80 +++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 59 insertions(+), 21 deletions(-) diff --git a/generators.py b/generators.py index f09aac7..e1282b7 100644 --- a/generators.py +++ b/generators.py @@ -3,6 +3,7 @@ from scipy.interpolate import CubicSpline from design import Design from geometry import rotate +from typing import Any def get_arc(R: float, start_deg: float, end_deg: float, n_points: int) -> np.ndarray: @@ -37,7 +38,11 @@ def generate_barriers(self) -> list[np.ndarray]: pass @abstractmethod - def generate_parameters(self) -> None: + def random_parameters(self) -> Any: + pass + + @abstractmethod + def set_parameters(self, params) -> None: pass @abstractmethod @@ -49,7 +54,8 @@ def name(self) -> str: return self.__class__.__name__ def random_barriers(self) -> list[np.ndarray]: - self.generate_parameters() + pars = self.random_parameters() + self.set_parameters(pars) return self.generate_barriers() def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: @@ -134,13 +140,20 @@ def _create_barrier( x_all, y_all = rotate(x_all, y_all, -45) return x_all, y_all - def generate_parameters(self) -> None: - # Parameterize the design - self.w_mins = self.w_mins_base + 0.5*np.random.random(4) + def random_parameters(self): + rand1 = 0.5*np.random.random(4) + rand2 = 2+3*np.random.random() + rand3 = 2+3*np.random.random() + rand4 = 2+2*np.random.random() + return rand1, rand2, rand3, rand4 + + def set_parameters(self, params) -> None: + rand1, rand2, rand3, rand4 = params + self.w_mins = self.w_mins_base + rand1 y_min0 = 13+4*np.random.random() - y_min1 = y_min0+self.w_mins[0]+2+3*np.random.random() - y_min2 = y_min1+self.w_mins[1]+2+3*np.random.random() - y_min3 = y_min2+self.w_mins[2]+2+2*np.random.random() + y_min1 = y_min0+self.w_mins[0] + rand2 + y_min2 = y_min1+self.w_mins[1] + rand3 + y_min3 = y_min2+self.w_mins[2] + rand4 self.w_mids = self.w_mins self.y_mins = np.array([y_min0, y_min1, y_min2, y_min3]) self.y_mids = self.y_mins + np.array([2, 1.5, 1, 0.5]) @@ -179,9 +192,15 @@ def __init__(self, design, r_stator_end, **kwargs): def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: pass - def generate_parameters(self) -> None: - self.phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) - self.phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) + def random_parameters(self): + phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) + phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) + return phis_inner, phis_outer + + def set_parameters(self, params) -> None: + phis_inner, phis_outer = params + self.phis_inner = phis_inner + self.phis_outer = phis_outer def generate_barriers(self) -> list[np.ndarray]: barriers = [] @@ -232,9 +251,15 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: ub = np.concatenate((self.phis_inner_max, [self.lam_max])) return lb, ub - def generate_parameters(self): - self.lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] - super().generate_parameters() + def random_parameters(self): + pars = super().random_parameters() + lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] + return *pars, lam + + def set_parameters(self, params) -> None: + phis_inner, phis_outer, lam = params + super().set_parameters((phis_inner, phis_outer)) + self.lam = lam def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y @@ -253,9 +278,15 @@ def __init__(self, design, r_stator_end, **kwargs): self.theta_max = 20 super().__init__(design, r_stator_end, **kwargs) - def generate_parameters(self): - self.theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] - super().generate_parameters() + def random_parameters(self): + pars = super().random_parameters() + theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] + return *pars, theta + + def set_parameters(self, params) -> None: + phis_inner, phis_outer, lam, theta = params + super().set_parameters((phis_inner, phis_outer, lam)) + self.theta = theta def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: x, y = super()._bezier_point(x, y, is_inner, order) @@ -280,10 +311,17 @@ def __init__(self, design, r_stator_end, **kwargs): self.lam_outer_max = 0.45 super().__init__(design, r_stator_end, **kwargs) - def generate_parameters(self): - self.lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] - self.lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] - super().generate_parameters() + def random_parameters(self): + pars = super().random_parameters() + lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] + lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] + return *pars, lam_inner, lam_outer + + def set_parameters(self, params) -> None: + phis_inner, phis_outer, lam_inner, lam_outer = params + super().set_parameters((phis_inner, phis_outer)) + self.lam_inner = lam_inner + self.lam_outer = lam_outer def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: lam = self.lam_inner if is_inner else self.lam_outer From 0590ed9ed61ae1bf42aaf846209cc1f051ef7aa9 Mon Sep 17 00:00:00 2001 From: Lukas Adam Date: Sun, 8 Mar 2026 14:19:52 +0100 Subject: [PATCH 28/55] Fixed bounds --- generators.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/generators.py b/generators.py index e1282b7..334a889 100644 --- a/generators.py +++ b/generators.py @@ -247,8 +247,8 @@ def __init__(self, design, r_stator_end, **kwargs): @property def bounds(self) -> tuple[np.ndarray, np.ndarray]: - lb = np.concatenate((self.phis_inner_min, [self.lam_min])) - ub = np.concatenate((self.phis_inner_max, [self.lam_max])) + lb = np.concatenate((self.phis_inner_min, self.phis_outer_min, [self.lam_min])) + ub = np.concatenate((self.phis_inner_max, self.phis_outer_max, [self.lam_max])) return lb, ub def random_parameters(self): From d06d94efb47292a4c98824adcba123f5eeaf82ad Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 15:43:41 +0000 Subject: [PATCH 29/55] Added bounds --- generators.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/generators.py b/generators.py index 334a889..cedf18b 100644 --- a/generators.py +++ b/generators.py @@ -140,6 +140,12 @@ def _create_barrier( x_all, y_all = rotate(x_all, y_all, -45) return x_all, y_all + @property + def bounds(self) -> tuple[np.ndarray, np.ndarray]: + lb = np.array([0, 0, 0, 0, 2, 2, 2]) + ub = np.array([0.5, 0.5, 0.5, 0.5, 5, 5, 4]) + return lb, ub + def random_parameters(self): rand1 = 0.5*np.random.random(4) rand2 = 2+3*np.random.random() From b3ed010177fe817d3f447614ea6cbaace869a8a3 Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 16:00:21 +0000 Subject: [PATCH 30/55] Removed save_barriers --- generators.py | 44 +++++--------------------------------------- 1 file changed, 5 insertions(+), 39 deletions(-) diff --git a/generators.py b/generators.py index cedf18b..422759e 100644 --- a/generators.py +++ b/generators.py @@ -1,6 +1,7 @@ from abc import ABC, abstractmethod import numpy as np from scipy.interpolate import CubicSpline +import pickle from design import Design from geometry import rotate from typing import Any @@ -26,6 +27,10 @@ def compute_r_max(design: Design, r_stator_end: float) -> float: airgap = design.mm_to_str("geom_params", "Airgap") return (dia_stator_gap / 2.0) - airgap - r_stator_end +def save_params(params, file_name: str) -> None: + with open(file_name, 'wb') as handle: + pickle.dump(params, handle, protocol=pickle.HIGHEST_PROTOCOL) + class BarrierGenerator(ABC): def __init__(self, offset: float, n_curve: int = 500, n_flat: int = 20) -> None: @@ -45,19 +50,10 @@ def random_parameters(self) -> Any: def set_parameters(self, params) -> None: pass - @abstractmethod - def save_barriers(self, file_name: str) -> None: - pass - @property def name(self) -> str: return self.__class__.__name__ - def random_barriers(self) -> list[np.ndarray]: - pars = self.random_parameters() - self.set_parameters(pars) - return self.generate_barriers() - def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: assert self.offset is not None check_barrier(barrier) @@ -173,15 +169,6 @@ def generate_barriers(self) -> list[np.ndarray]: xy_all = np.vstack((x_all, y_all)).T barriers.append(xy_all) return barriers - - def save_barriers(self, file_name: str) -> None: - np.savez(file_name, - y_mins=self.y_mins, - w_mins=self.w_mins, - y_mids=self.y_mids, - w_mids=self.w_mids, - thetas=self.thetas, - w_maxs=self.w_maxs) class AbstractHacklGenerator(BarrierGenerator): @@ -270,13 +257,6 @@ def set_parameters(self, params) -> None: def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y - def save_barriers(self, file_name: str): - np.savez(file_name, - phis_inner=self.phis_inner, - phis_outer=self.phis_outer, - lam=self.lam) - - class HacklGenerator_OneLambdaTheta(HacklGenerator_OneLambda): def __init__(self, design, r_stator_end, **kwargs): @@ -301,13 +281,6 @@ def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple return xx[0], yy[0] return x, y - def save_barriers(self, file_name: str): - np.savez(file_name, - phis_inner=self.phis_inner, - phis_outer=self.phis_outer, - lam=self.lam, - theta=self.theta) - class HacklGenerator_TwoLambdas(AbstractHacklGenerator): def __init__(self, design, r_stator_end, **kwargs): @@ -332,10 +305,3 @@ def set_parameters(self, params) -> None: def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: lam = self.lam_inner if is_inner else self.lam_outer return lam * x + (1 - lam) * y, y - - def save_barriers(self, file_name: str): - np.savez(file_name, - phis_inner=self.phis_inner, - phis_outer=self.phis_outer, - lam_inner=self.lam_inner, - lam_outer=self.lam_outer) \ No newline at end of file From eccfacfa99399a7afcdbf858544c003b4fff727c Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 16:01:07 +0000 Subject: [PATCH 31/55] Removed save_barriers --- run.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/run.py b/run.py index c89b15b..48d7c81 100644 --- a/run.py +++ b/run.py @@ -4,7 +4,7 @@ import numpy as np import pandas as pd from design import Design -from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta +from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta, save_params from geometry import analyze_results, plot_barriers project_name = "SynRM_test" @@ -51,7 +51,9 @@ metadata = pd.DataFrame() for i in range(0, n_designs): for generator in generators: - barriers = generator.random_barriers() + params = generator.random_parameters() + generator.set_parameters(params) + barriers = generator.generate_barriers() barriers = generator.split_barriers(barriers) design.add_rotor() @@ -74,7 +76,7 @@ title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %' file_name = f'{path_results}/design_{generator.name}_{i}' plot_barriers(barriers, design, title=title, file_name = f"{file_name}.png") - generator.save_barriers(f"{file_name}.npz") + save_params(params, f"{file_name}.pkl") metadata_new = { 'method': generator.name, From 9dd950e1f146a4b516a95e83c01118d0db58590d Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 16:10:08 +0000 Subject: [PATCH 32/55] Added X_to_params --- generators.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/generators.py b/generators.py index 422759e..85a677d 100644 --- a/generators.py +++ b/generators.py @@ -50,6 +50,10 @@ def random_parameters(self) -> Any: def set_parameters(self, params) -> None: pass + @abstractmethod + def X_to_params(self, X: np.ndarray) -> tuple[Any, ...]: + pass + @property def name(self) -> str: return self.__class__.__name__ @@ -162,6 +166,9 @@ def set_parameters(self, params) -> None: self.thetas = [1, 8, 15, 21] self.w_maxs = self.w_mins - np.array([0.5, 0.5, 0.5, 0]) + def X_to_params(self, X: np.ndarray): + return X[:4], X[4], X[5], X[6] + def generate_barriers(self) -> list[np.ndarray]: barriers = [] for args in zip(self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs): @@ -254,6 +261,9 @@ def set_parameters(self, params) -> None: super().set_parameters((phis_inner, phis_outer)) self.lam = lam + def X_to_params(self, X: np.ndarray): + return X[:3], X[:6], X[6] + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y @@ -274,6 +284,9 @@ def set_parameters(self, params) -> None: super().set_parameters((phis_inner, phis_outer, lam)) self.theta = theta + def X_to_params(self, X: np.ndarray): + return X[:3], X[:6], X[6], X[7] + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: x, y = super()._bezier_point(x, y, is_inner, order) if order < 2: @@ -302,6 +315,9 @@ def set_parameters(self, params) -> None: self.lam_inner = lam_inner self.lam_outer = lam_outer + def X_to_params(self, X: np.ndarray): + return X[:3], X[:6], X[6], X[7] + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: lam = self.lam_inner if is_inner else self.lam_outer return lam * x + (1 - lam) * y, y From 52a3088ae9d8a4c69c72d20d8cdf8f3b8301b7fc Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 16:10:45 +0000 Subject: [PATCH 33/55] Gitignore updated --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 218a65e..6caafd6 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ data/ results*/ nohup.out +.vscode/settings.json # Byte-compiled / optimized / DLL files __pycache__/ From 8e008a24812b90bfb8f3705b0cf0701f10f28d5d Mon Sep 17 00:00:00 2001 From: sadda Date: Sun, 8 Mar 2026 17:08:07 +0000 Subject: [PATCH 34/55] Fixed X_to_params --- generators.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/generators.py b/generators.py index 85a677d..a426fd9 100644 --- a/generators.py +++ b/generators.py @@ -262,7 +262,7 @@ def set_parameters(self, params) -> None: self.lam = lam def X_to_params(self, X: np.ndarray): - return X[:3], X[:6], X[6] + return X[:3], X[3:6], X[6] def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y @@ -285,7 +285,7 @@ def set_parameters(self, params) -> None: self.theta = theta def X_to_params(self, X: np.ndarray): - return X[:3], X[:6], X[6], X[7] + return X[:3], X[3:6], X[6], X[7] def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: x, y = super()._bezier_point(x, y, is_inner, order) @@ -316,7 +316,7 @@ def set_parameters(self, params) -> None: self.lam_outer = lam_outer def X_to_params(self, X: np.ndarray): - return X[:3], X[:6], X[6], X[7] + return X[:3], X[3:6], X[6], X[7] def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: lam = self.lam_inner if is_inner else self.lam_outer From 1e653806a617029a2184c61476b3fcf00174e61b Mon Sep 17 00:00:00 2001 From: sadda Date: Mon, 9 Mar 2026 05:07:49 +0000 Subject: [PATCH 35/55] Fixed random param --- generators.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/generators.py b/generators.py index a426fd9..1c1662f 100644 --- a/generators.py +++ b/generators.py @@ -142,24 +142,25 @@ def _create_barrier( @property def bounds(self) -> tuple[np.ndarray, np.ndarray]: - lb = np.array([0, 0, 0, 0, 2, 2, 2]) - ub = np.array([0.5, 0.5, 0.5, 0.5, 5, 5, 4]) + lb = np.array([0, 0, 0, 0, 13, 2, 2, 2]) + ub = np.array([0.5, 0.5, 0.5, 0.5, 17, 5, 5, 4]) return lb, ub def random_parameters(self): rand1 = 0.5*np.random.random(4) - rand2 = 2+3*np.random.random() + rand2 = 13+4*np.random.random() rand3 = 2+3*np.random.random() - rand4 = 2+2*np.random.random() - return rand1, rand2, rand3, rand4 + rand4 = 2+3*np.random.random() + rand5 = 2+2*np.random.random() + return rand1, rand2, rand3, rand4, rand5 def set_parameters(self, params) -> None: - rand1, rand2, rand3, rand4 = params + rand1, rand2, rand3, rand4, rand5 = params self.w_mins = self.w_mins_base + rand1 - y_min0 = 13+4*np.random.random() - y_min1 = y_min0+self.w_mins[0] + rand2 - y_min2 = y_min1+self.w_mins[1] + rand3 - y_min3 = y_min2+self.w_mins[2] + rand4 + y_min0 = rand2 + y_min1 = y_min0+self.w_mins[0] + rand3 + y_min2 = y_min1+self.w_mins[1] + rand4 + y_min3 = y_min2+self.w_mins[2] + rand5 self.w_mids = self.w_mins self.y_mins = np.array([y_min0, y_min1, y_min2, y_min3]) self.y_mids = self.y_mins + np.array([2, 1.5, 1, 0.5]) @@ -167,7 +168,7 @@ def set_parameters(self, params) -> None: self.w_maxs = self.w_mins - np.array([0.5, 0.5, 0.5, 0]) def X_to_params(self, X: np.ndarray): - return X[:4], X[4], X[5], X[6] + return X[:4], X[4], X[5], X[6], X[7] def generate_barriers(self) -> list[np.ndarray]: barriers = [] From d6a0fec15407acad42efdfa4dfdb741fb9f47cc3 Mon Sep 17 00:00:00 2001 From: sadda Date: Mon, 9 Mar 2026 05:08:06 +0000 Subject: [PATCH 36/55] Changed iterations --- run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/run.py b/run.py index 48d7c81..0c4ef63 100644 --- a/run.py +++ b/run.py @@ -15,7 +15,7 @@ os.makedirs(path, exist_ok=True) file_name_aedt = f'{path_data}/{project_name}.aedt' plot_design = True -n_designs = 100 +n_designs = 20 # Define constants AEDT_VERSION = "2024.1" From de13ca8eb4aa5219061b574eccf58c1767329214 Mon Sep 17 00:00:00 2001 From: sadda Date: Mon, 9 Mar 2026 18:30:43 +0000 Subject: [PATCH 37/55] Added script --- run_optimization.py | 167 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 run_optimization.py diff --git a/run_optimization.py b/run_optimization.py new file mode 100644 index 0000000..e98da44 --- /dev/null +++ b/run_optimization.py @@ -0,0 +1,167 @@ +# export ANSYSEM_ROOT241=/data/AnsysEM/v241/Linux64 + +from collections.abc import Iterable +import os +import numpy as np +import pandas as pd +import pickle +from design import Design +from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta +from geometry import analyze_results, plot_barriers +import torch +torch.set_default_dtype(torch.float64) +from torch import Tensor +from botorch.utils.transforms import normalize, unnormalize +import torch + +from botorch.models import SingleTaskGP +from botorch import fit_gpytorch_mll +from gpytorch.mlls import ExactMarginalLogLikelihood + +from botorch.acquisition.multi_objective.logei import qLogExpectedHypervolumeImprovement +from botorch.utils.multi_objective.box_decompositions import NondominatedPartitioning +from botorch.utils.multi_objective.pareto import is_non_dominated +from botorch.optim import optimize_acqf + +def objective_transform(TorAvg, TorRippleRms): + if pd.isnull(TorAvg): + return -99999, -99999 + else: + return TorAvg, -TorRippleRms/100 + +def objective(Xs: Tensor, *args) -> Tensor: + vals = [objective_single(X, *args) for X in Xs] + return torch.stack(vals, dim=0) + +def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: + X = unnormalize(X, bounds) + params = generator.X_to_params(X.numpy()) + + generator.set_parameters(params) + barriers = generator.generate_barriers() + barriers = generator.split_barriers(barriers) + + design.add_rotor() + + for barrier in barriers: + design.add_rotor_barrier(barrier) + + # Compute the torque + Tor = design.compute(NUM_CORES) + if Tor is None: + TorAvg, TorRippleRms = np.nan, np.nan + else: + TorAvg, _, TorRippleRms = analyze_results(Tor) + + # Delete the rotor + design.delete_rotor() + + f1, f2 = objective_transform(TorAvg, TorRippleRms) + return torch.tensor([f1, f2]) + +def init_points(root, method): + results = pd.read_csv(f"{root}/metadata.csv") + results = results[~results["T"].isnull()] + results = results[results["method"] == method] + results["path"] = [f"results/design_{row['method']}_{row['design']}.pkl" for _, row in results.iterrows()] + + Xs, Ys = [], [] + for _, r in results.iterrows(): + with open(r["path"], 'rb') as f: + params = pickle.load(f) + + X = [] + for x in params: + if isinstance(x, Iterable): + for y in x: + X.append(y) + else: + X.append(x) + Y = objective_transform(r["T"], r["ripple"]) + + Xs.append(torch.Tensor(X)) + Ys.append(torch.Tensor(Y)) + + return torch.stack(Xs), torch.stack(Ys) + +project_name = "SynRM_test" +design_name = "Design01" +path_data = os.path.join(os.getcwd(), 'data') +os.makedirs(path_data, exist_ok=True) +file_name_aedt = f'{path_data}/{project_name}.aedt' + +# Define constants +AEDT_VERSION = "2024.1" +NUM_CORES = 4 +NG_MODE = True #non-graphical mode +CLS_EXIT = True #close on exit + +if not os.path.exists(file_name_aedt): + design = Design.create( + project_name, design_name, file_name_aedt, + version=AEDT_VERSION, + non_graphical=NG_MODE, + new_desktop=False, + close_on_exit=CLS_EXIT, + ) +else: + design = Design.load( + file_name_aedt, + version=AEDT_VERSION, + non_graphical=NG_MODE, + new_desktop=False, + close_on_exit=CLS_EXIT, + ) + +n_iters = 40 +batch_size = 4 +r_stator_end = 0.7 +offset = 0.7 / 2 +# generator = HacklGenerator_OneLambda(design, r_stator_end, offset=offset) +generator = FourStupid(design, r_stator_end, offset=offset) +bounds = torch.from_numpy(np.vstack(generator.bounds)) + +root_init = "results" +method = generator.__class__.__name__ +train_X, train_Y = init_points(root_init, method) +train_X = normalize(train_X, bounds) +bounds_normalized = normalize(bounds, bounds) +ref_point = torch.tensor([3.8,-0.3]) + +def objective_lambda(Xs): + return objective(Xs, design, generator, bounds, NUM_CORES) + +for _ in range(n_iters): + # Fit surrogate + model = SingleTaskGP(train_X, train_Y) + mll = ExactMarginalLogLikelihood(model.likelihood, model) + fit_gpytorch_mll(mll) + + # Compute Pareto front + pareto_Y = train_Y[is_non_dominated(train_Y)] + partitioning = NondominatedPartitioning( + ref_point=ref_point, + Y=pareto_Y + ) + + # Optimize acquisition function to select candidate points + acq = qLogExpectedHypervolumeImprovement( + model=model, + ref_point=ref_point.tolist(), + partitioning=partitioning, + ) + candidates, _ = optimize_acqf( + acq_function=acq, + bounds=bounds_normalized, + q=batch_size, + num_restarts=10, + raw_samples=128, + ) + + # Evaluate candidates + new_Y = objective_lambda(candidates) + train_X = torch.cat([train_X, candidates]) + train_Y = torch.cat([train_Y, new_Y]) + np.savez(f"results_{method}.npz", train_X=unnormalize(train_X, bounds), train_Y=train_Y) + +design.close_project() \ No newline at end of file From 9a22ad620158fe30ae0a11d019395dabb9e07ae4 Mon Sep 17 00:00:00 2001 From: sadda Date: Mon, 9 Mar 2026 18:31:00 +0000 Subject: [PATCH 38/55] Updated gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 6caafd6..687c547 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ data/ results*/ +results_*.npz nohup.out .vscode/settings.json From 05e76ce1fcf878804bee5125045839fae766d284 Mon Sep 17 00:00:00 2001 From: sadda Date: Mon, 9 Mar 2026 18:32:51 +0000 Subject: [PATCH 39/55] Ruff --- design.py | 165 +++++++++++++++++++++++++++----------------- generators.py | 153 ++++++++++++++++++++++++---------------- geometry.py | 20 +++--- run.py | 44 +++++++----- run_optimization.py | 52 +++++++++----- 5 files changed, 267 insertions(+), 167 deletions(-) diff --git a/design.py b/design.py index 819e974..9cae6b8 100644 --- a/design.py +++ b/design.py @@ -1,20 +1,22 @@ -import os import numpy as np from ansys.aedt.core import Desktop, Maxwell2d from ansys.aedt.core.modeler.modeler_2d import Modeler2D + class Design: def __init__(self, m2d: Maxwell2d) -> None: self.set_parameters() self.m2d = m2d @classmethod - def create(cls, project_name: str, design_name: str, file_name: str, **kwargs) -> "Design": + def create( + cls, project_name: str, design_name: str, file_name: str, **kwargs + ) -> "Design": m2d = Maxwell2d( project=project_name, design=design_name, solution_type="TransientXY", - **kwargs + **kwargs, ) obj = cls(m2d) obj.create_stator() @@ -30,9 +32,9 @@ def load(cls, file_name: str, **kwargs) -> "Design": m2d.set_active_design("Design01") return cls(m2d) - + def set_parameters(self) -> None: - #materials + # materials self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" # main definitions @@ -43,7 +45,7 @@ def set_parameters(self) -> None: "SlotNumber": "36", "SlotType": "3", "DiaShaft": "25mm", - "StackLength": "85mm" + "StackLength": "85mm", } # stator slot self.slot_params = { @@ -63,40 +65,57 @@ def set_parameters(self) -> None: "CoilPitch": "9", # coil pitch in slots "SlotLiner": "0.3mm", "SpaceLayers": "0.2mm", - "Nc": "68" # turns per coil + "Nc": "68", # turns per coil } # model parameters PolePairs = 2 - f = 50 #[Hz] - RotSpeed = 60*f/PolePairs #[rpm] + f = 50 # [Hz] + RotSpeed = 60 * f / PolePairs # [rpm] self.mod_params = { "Poles": f"2*{PolePairs}", "ModelLength": "85mm", "SymmetryFactor": "Poles", "StatorSkewAngle": "0deg", } - #operation parameters + # operation parameters self.oper_params = { "Im": "1.5*sqrt(2)A", - "epsI": "pi/4", #current angle + "epsI": "pi/4", # current angle "InitPos": "-30deg", "f": f"{f}Hz", "RotSpeed": f"{RotSpeed}rpm", - "Nper": "1/6" - "", #number of included periods - "PointPer": "101" #number of time points per period + "Nper": "1/6", # number of included periods + "PointPer": "101", # number of time points per period } self.rot_points = [ - ["DiaShaft/2*cos(360deg/SymmetryFactor)", "DiaShaft/2*sin(360deg/SymmetryFactor)", "0mm"], - ["DiaShaft/2*cos(360deg/(2*SymmetryFactor))", "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", "0mm"], + [ + "DiaShaft/2*cos(360deg/SymmetryFactor)", + "DiaShaft/2*sin(360deg/SymmetryFactor)", + "0mm", + ], + [ + "DiaShaft/2*cos(360deg/(2*SymmetryFactor))", + "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", + "0mm", + ], ["DiaShaft/2", "0mm", "0mm"], ["DiaStatorGap/2-Airgap", "0mm", "0mm"], - ["(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", "0mm"], - ["(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", "0mm"], + [ + "(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", + "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", + "0mm", + ], + [ + "(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", + "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", + "0mm", + ], ] self.setup_name = "Setup1" - self.rotor_r_min = self.mm_to_str('geom_params', 'DiaShaft')/2 - self.rotor_r_max = self.mm_to_str('geom_params', 'DiaStatorGap')/2 - self.mm_to_str('geom_params', 'Airgap') + self.rotor_r_min = self.mm_to_str("geom_params", "DiaShaft") / 2 + self.rotor_r_max = self.mm_to_str( + "geom_params", "DiaStatorGap" + ) / 2 - self.mm_to_str("geom_params", "Airgap") def create_stator(self) -> None: m2d = self.m2d @@ -141,7 +160,7 @@ def create_stator(self) -> None: name="Shaft", ) - #motion setup + # motion setup m2d.assign_rotate_motion( assignment="Band", coordinate_system="Global", @@ -199,10 +218,14 @@ def create_stator(self) -> None: stator_id.color = (192, 192, 192) # rgb stator_id.transparency = 0.0 - #Winding + # Winding coil_id = modeler.create_rectangle( origin=["DiaStatorGap/2+Hs0+Hs1+SlotLiner", "-(Bs1/2-SlotLiner)", 0], - sizes=["(Hs2+Rs/2-2*SlotLiner)/Layers-(Layers-1)*SpaceLayers/2", "Bs1-2*SlotLiner", 0], + sizes=[ + "(Hs2+Rs/2-2*SlotLiner)/Layers-(Layers-1)*SpaceLayers/2", + "Bs1-2*SlotLiner", + 0, + ], name="Coil", material="copper", ) @@ -210,7 +233,10 @@ def create_stator(self) -> None: coil_id.transparency = 0.0 modeler.rotate(assignment=coil_id, axis="Z", angle="360deg/SlotNumber/2") coil_id.duplicate_around_axis( - axis="Z", angle="360deg/SlotNumber", clones="CoilPitch", create_new_objects=True + axis="Z", + angle="360deg/SlotNumber", + clones="CoilPitch", + create_new_objects=True, ) id_coils = modeler.get_objects_w_string(string_name="Coil", case_sensitive=True) @@ -260,16 +286,14 @@ def create_stator(self) -> None: ], assignment="Region", ) - m2d.assign_vector_potential( - assignment=id_bc_az, vector_value=0, boundary="A0" - ) + m2d.assign_vector_potential(assignment=id_bc_az, vector_value=0, boundary="A0") # Excitations I_A = "Im * cos(2*pi*f*time+epsI)" I_B = "Im * cos(2*pi*f*time-120deg+epsI)" I_C = "Im * cos(2*pi*f*time-240deg+epsI)" - #Define phase windings + # Define phase windings m2d.assign_coil( assignment=["Coil"], conductors_number="Nc", @@ -296,9 +320,7 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseA", ) - m2d.add_winding_coils( - assignment="PhaseA", coils=["CS1", "CS2", "CS3"] - ) + m2d.add_winding_coils(assignment="PhaseA", coils=["CS1", "CS2", "CS3"]) m2d.assign_coil( assignment=["Coil_6"], conductors_number="Nc", @@ -325,9 +347,7 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseB", ) - m2d.add_winding_coils( - assignment="PhaseB", coils=["CS7", "CS8", "CS9"] - ) + m2d.add_winding_coils(assignment="PhaseB", coils=["CS7", "CS8", "CS9"]) m2d.assign_coil( assignment=["Coil_3"], conductors_number="Nc", @@ -354,9 +374,7 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseC", ) - m2d.add_winding_coils( - assignment="PhaseC", coils=["CS4", "CS5", "CS6"] - ) + m2d.add_winding_coils(assignment="PhaseC", coils=["CS4", "CS5", "CS6"]) # Mesh operation m2d.mesh.assign_length_mesh( @@ -374,18 +392,18 @@ def create_stator(self) -> None: name="stator", ) - #core loss + # core loss m2d.set_core_losses("Stator", core_loss_on_field=False) - #inductance calculation + # inductance calculation m2d.change_inductance_computation( compute_transient_inductance=True, incremental_matrix=False ) - #model depth + # model depth m2d.model_depth = "StackLength" - #symmetry + # symmetry m2d.change_symmetry_multiplier("SymmetryFactor") - # Calculation setup + # Calculation setup setup = m2d.create_setup(name=self.setup_name) setup.props["StopTime"] = "Nper/f" setup.props["TimeStep"] = "1/(f*(PointPer-1))" @@ -396,7 +414,7 @@ def create_stator(self) -> None: setup.update() m2d.validate_simple() - #ooutput variables + # ooutput variables output_vars = { "pos": "(Moving1.Position -InitPos) * Poles/2", "cos0": "cos(pos)", @@ -424,21 +442,34 @@ def create_stator(self) -> None: for k, v in output_vars.items(): m2d.create_output_variable(k, v) - #Definitions for plots + # Definitions for plots post_params = { # reports - ("InducedVoltage(PhaseA)","InducedVoltage(PhaseB)","InducedVoltage(PhaseC)"): "InducedVoltage", + ( + "InducedVoltage(PhaseA)", + "InducedVoltage(PhaseB)", + "InducedVoltage(PhaseC)", + ): "InducedVoltage", ("Moving1.Torque"): "Torque", - ("InputCurrent(PhaseA)","InputCurrent(PhaseB)","InputCurrent(PhaseC)"): "Current", - ("FluxLinkage(PhaseA)","FluxLinkage(PhaseB)","FluxLinkage(PhaseC)",): "FluxLinkage", + ( + "InputCurrent(PhaseA)", + "InputCurrent(PhaseB)", + "InputCurrent(PhaseC)", + ): "Current", + ( + "FluxLinkage(PhaseA)", + "FluxLinkage(PhaseB)", + "FluxLinkage(PhaseC)", + ): "FluxLinkage", ("I_d", "I_q"): "Current_dq", ("Flux_d", "Flux_q"): "FluxLinkage_dq", ("Ui_d", "Ui_q"): "InducedVoltage_dq", ("L_d", "L_q"): "Inductance_dq", - } - #Create Report + # Create Report for k, v in post_params.items(): - expressions = list(k) if isinstance(k, tuple) else [k] #if multiple report, use list(k). Else, use k + expressions = ( + list(k) if isinstance(k, tuple) else [k] + ) # if multiple report, use list(k). Else, use k report = m2d.post.create_report( expressions=expressions, setup_sweep_name="", @@ -459,7 +490,10 @@ def add_rotor(self) -> None: assert isinstance(modeler, Modeler2D) rotor_id = modeler.create_polyline( - points=self.rot_points, segment_type=[ "Arc","Line", "Arc"], cover_surface=True, name="Rotor" + points=self.rot_points, + segment_type=["Arc", "Line", "Arc"], + cover_surface=True, + name="Rotor", ) self.rotor_id = rotor_id rotor_id.material_name = self.Fe @@ -473,14 +507,19 @@ def add_rotor_barrier(self, barrier_points, segment_type=None) -> None: # Round it for Ansys barrier_points = np.round(barrier_points, 6) - # Potentially add the z axis + # Potentially add the z axis if barrier_points.shape[1] == 2: - barrier_points = np.hstack((barrier_points, np.zeros((len(barrier_points),1)))) - + barrier_points = np.hstack( + (barrier_points, np.zeros((len(barrier_points), 1))) + ) + # Convert them into a string format and interpolate points_str = [[str(y) for y in x] for x in barrier_points] barrier_id = modeler.create_polyline( - points=points_str, segment_type=segment_type, cover_surface=True, name="Barrier" + points=points_str, + segment_type=segment_type, + cover_surface=True, + name="Barrier", ) # Remove the barrier @@ -499,7 +538,7 @@ def compute(self, NUM_CORES: int = 1): maximum_elements=None, name="rotor", ) - #core loss rotor + # core loss rotor m2d.set_core_losses("Rotor", core_loss_on_field=False) # Analyze @@ -512,14 +551,14 @@ def compute(self, NUM_CORES: int = 1): return solutions.data_magnitude() except AttributeError: return None - + def delete_rotor(self) -> None: assert isinstance(self.m2d.modeler, Modeler2D) self.m2d.modeler.delete(self.rotor_id) def save_design(self, file_name: str, **kwargs) -> None: - show = kwargs.pop('show', False) - view = kwargs.pop('view', 'xy') + show = kwargs.pop("show", False) + view = kwargs.pop("view", "xy") self.m2d.plot(show=show, output_file=file_name, view=view) def save_project(self, file_name: str | None = None) -> None: @@ -527,12 +566,12 @@ def save_project(self, file_name: str | None = None) -> None: self.m2d.save_project() else: self.m2d.save_project(file_name) - + def close_project(self) -> None: self.m2d.close_desktop() - + def mm_to_str(self, var, field) -> float: val = getattr(self, var)[field] - if not val.endswith('mm'): - raise Exception('val must end with mm') + if not val.endswith("mm"): + raise Exception("val must end with mm") return float(val[:-2]) diff --git a/generators.py b/generators.py index 1c1662f..e1e5f32 100644 --- a/generators.py +++ b/generators.py @@ -15,20 +15,24 @@ def get_arc(R: float, start_deg: float, end_deg: float, n_points: int) -> np.nda return np.column_stack((x, y)) + def signed_distance(P: np.ndarray, a: float, b: float, c: float) -> np.ndarray: - return a * P[:,0] + b * P[:,1] + c + return a * P[:, 0] + b * P[:, 1] + c + def check_barrier(barrier: np.ndarray) -> None: if not np.array_equal(barrier[0], barrier[-1]): raise ValueError("Barrier must start and end at the same point") + def compute_r_max(design: Design, r_stator_end: float) -> float: dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") airgap = design.mm_to_str("geom_params", "Airgap") return (dia_stator_gap / 2.0) - airgap - r_stator_end + def save_params(params, file_name: str) -> None: - with open(file_name, 'wb') as handle: + with open(file_name, "wb") as handle: pickle.dump(params, handle, protocol=pickle.HIGHEST_PROTOCOL) @@ -37,7 +41,7 @@ def __init__(self, offset: float, n_curve: int = 500, n_flat: int = 20) -> None: self.offset = offset self.n_curve = n_curve self.n_flat = n_flat - + @abstractmethod def generate_barriers(self) -> list[np.ndarray]: pass @@ -64,7 +68,7 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: lam = np.linspace(0, 1, self.n_flat)[:, None] barriers = [] - for a, b, c in zip([1., -1.], [-1., 1.], [-self.offset, -self.offset]): + for a, b, c in zip([1.0, -1.0], [-1.0, 1.0], [-self.offset, -self.offset]): f = signed_distance(barrier, a, b, c) f_is_positive = f > 0 mask_positive = np.where(f_is_positive)[0] @@ -72,7 +76,7 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: barrier_new = [barrier[mask_positive[0]]] for i in range(len(mask_positive) - 1): i1 = mask_positive[i] - i2 = mask_positive[i+1] + i2 = mask_positive[i + 1] if i2 - i1 == 1: barrier_new.append(barrier[i2]) @@ -92,8 +96,11 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: return [x for barrier in barriers for x in self.split_barrier(barrier)] + class FourStupid(BarrierGenerator): - def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwargs) -> None: + def __init__( + self, design, r_stator_end, der1=1.0, der2=1.0, symmetric=True, **kwargs + ) -> None: self.R = compute_r_max(design, r_stator_end) self.design = design self.der1 = der1 @@ -103,31 +110,31 @@ def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwa super().__init__(**kwargs) def _create_barrier( - self, - y_min, - w_min, - y_mid, - w_mid, - theta, - w_max, - ): - - theta1 = (theta+45) / 180*np.pi - x_max1 = self.R*np.cos(theta1) - y_max1 = self.R*np.sin(theta1) - - x1 = [0, x_max1/2, x_max1] + self, + y_min, + w_min, + y_mid, + w_mid, + theta, + w_max, + ): + + theta1 = (theta + 45) / 180 * np.pi + x_max1 = self.R * np.cos(theta1) + y_max1 = self.R * np.sin(theta1) + + x1 = [0, x_max1 / 2, x_max1] y1 = [y_min, y_mid, y_max1] - f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,self.der1))) + f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1, self.der1))) theta2 = theta1 + w_max / self.design.rotor_r_max - x_max2 = self.R*np.cos(theta2) - y_max2 = self.R*np.sin(theta2) + x_max2 = self.R * np.cos(theta2) + y_max2 = self.R * np.sin(theta2) - s = w_mid / np.sqrt(1 + self.der1**2/4) - x2 = [0, x1[1]-s*self.der1/2, x_max2] - y2 = [y_min+w_min, y1[1]+s, y_max2] - f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,self.der2))) + s = w_mid / np.sqrt(1 + self.der1**2 / 4) + x2 = [0, x1[1] - s * self.der1 / 2, x_max2] + y2 = [y_min + w_min, y1[1] + s, y_max2] + f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1, self.der2))) x_interp1 = np.linspace(x1[0], x1[-1], self.n_curve) x_interp2 = np.linspace(x2[0], x2[-1], self.n_curve) @@ -147,20 +154,20 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: return lb, ub def random_parameters(self): - rand1 = 0.5*np.random.random(4) - rand2 = 13+4*np.random.random() - rand3 = 2+3*np.random.random() - rand4 = 2+3*np.random.random() - rand5 = 2+2*np.random.random() + rand1 = 0.5 * np.random.random(4) + rand2 = 13 + 4 * np.random.random() + rand3 = 2 + 3 * np.random.random() + rand4 = 2 + 3 * np.random.random() + rand5 = 2 + 2 * np.random.random() return rand1, rand2, rand3, rand4, rand5 def set_parameters(self, params) -> None: rand1, rand2, rand3, rand4, rand5 = params self.w_mins = self.w_mins_base + rand1 y_min0 = rand2 - y_min1 = y_min0+self.w_mins[0] + rand3 - y_min2 = y_min1+self.w_mins[1] + rand4 - y_min3 = y_min2+self.w_mins[2] + rand5 + y_min1 = y_min0 + self.w_mins[0] + rand3 + y_min2 = y_min1 + self.w_mins[1] + rand4 + y_min3 = y_min2 + self.w_mins[2] + rand5 self.w_mids = self.w_mins self.y_mins = np.array([y_min0, y_min1, y_min2, y_min3]) self.y_mids = self.y_mins + np.array([2, 1.5, 1, 0.5]) @@ -172,30 +179,38 @@ def X_to_params(self, X: np.ndarray): def generate_barriers(self) -> list[np.ndarray]: barriers = [] - for args in zip(self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs): + for args in zip( + self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs + ): x_all, y_all = self._create_barrier(*args) xy_all = np.vstack((x_all, y_all)).T barriers.append(xy_all) return barriers - + class AbstractHacklGenerator(BarrierGenerator): def __init__(self, design, r_stator_end, **kwargs): self.R = compute_r_max(design, r_stator_end) - self.phis_inner_min = np.asarray([4., 17.3, 30.6]) + self.phis_inner_min = np.asarray([4.0, 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) - self.phis_outer_min = np.asarray([10.6, 24., 37.3]) - self.phis_outer_max = np.asarray([16.8, 30.2, 44.]) + self.phis_outer_min = np.asarray([10.6, 24.0, 37.3]) + self.phis_outer_max = np.asarray([16.8, 30.2, 44.0]) self.n_barriers = len(self.phis_inner_min) super().__init__(**kwargs) @abstractmethod - def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: + def _bezier_point( + self, x: float, y: float, is_inner: bool, order: int + ) -> tuple[float, float]: pass def random_parameters(self): - phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) - phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) + phis_inner = self.phis_inner_min + ( + self.phis_inner_max - self.phis_inner_min + ) * np.random.rand(self.n_barriers) + phis_outer = self.phis_outer_min + ( + self.phis_outer_max - self.phis_outer_min + ) * np.random.rand(self.n_barriers) return phis_inner, phis_outer def set_parameters(self, params) -> None: @@ -205,17 +220,21 @@ def set_parameters(self, params) -> None: def generate_barriers(self) -> list[np.ndarray]: barriers = [] - for i, (phi_inner, phi_outer) in enumerate(zip(self.phis_inner, self.phis_outer)): + for i, (phi_inner, phi_outer) in enumerate( + zip(self.phis_inner, self.phis_outer) + ): # Generate the longer part of the barrier pts_outer = self._get_bezier_curve(phi_outer, False, i) pts_inner = self._get_bezier_curve(phi_inner, True, i) - + # Generate the arcs at the end arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner, self.n_flat) arc_bottom = get_arc(self.R, phi_inner, phi_outer, self.n_flat) - + # Merge them together - barrier = np.concatenate((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) + barrier = np.concatenate( + (pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:]) + ) barriers.append(barrier) return barriers @@ -233,8 +252,8 @@ def _get_bezier_curve(self, phi_deg, is_inner, i): # Formula 33: Polynomial expansion to generate continuous points z_vals = np.linspace(0, 1, self.n_curve)[:, None] return ( - (1 - z_vals)**3 * r0 - + 3 * (1 - z_vals)**2 * z_vals * r1 + (1 - z_vals) ** 3 * r0 + + 3 * (1 - z_vals) ** 2 * z_vals * r1 + 3 * (1 - z_vals) * z_vals**2 * r2 + z_vals**3 * r3 ) @@ -251,12 +270,12 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: lb = np.concatenate((self.phis_inner_min, self.phis_outer_min, [self.lam_min])) ub = np.concatenate((self.phis_inner_max, self.phis_outer_max, [self.lam_max])) return lb, ub - + def random_parameters(self): pars = super().random_parameters() - lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] + lam = self.lam_min + (self.lam_max - self.lam_min) * np.random.rand(1)[0] return *pars, lam - + def set_parameters(self, params) -> None: phis_inner, phis_outer, lam = params super().set_parameters((phis_inner, phis_outer)) @@ -264,8 +283,10 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6] - - def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: + + def _bezier_point( + self, x: float, y: float, is_inner: bool, order: int + ) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y @@ -277,7 +298,9 @@ def __init__(self, design, r_stator_end, **kwargs): def random_parameters(self): pars = super().random_parameters() - theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] + theta = ( + self.theta_min + (self.theta_max - self.theta_min) * np.random.rand(1)[0] + ) return *pars, theta def set_parameters(self, params) -> None: @@ -287,8 +310,10 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6], X[7] - - def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: + + def _bezier_point( + self, x: float, y: float, is_inner: bool, order: int + ) -> tuple[float, float]: x, y = super()._bezier_point(x, y, is_inner, order) if order < 2: xx, yy = rotate(x, y, self.theta) @@ -306,8 +331,14 @@ def __init__(self, design, r_stator_end, **kwargs): def random_parameters(self): pars = super().random_parameters() - lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] - lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] + lam_inner = ( + self.lam_inner_min + + (self.lam_inner_max - self.lam_inner_min) * np.random.rand(1)[0] + ) + lam_outer = ( + self.lam_outer_min + + (self.lam_outer_max - self.lam_outer_min) * np.random.rand(1)[0] + ) return *pars, lam_inner, lam_outer def set_parameters(self, params) -> None: @@ -319,6 +350,8 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6], X[7] - def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: + def _bezier_point( + self, x: float, y: float, is_inner: bool, order: int + ) -> tuple[float, float]: lam = self.lam_inner if is_inner else self.lam_outer return lam * x + (1 - lam) * y, y diff --git a/geometry.py b/geometry.py index 328ba07..3e659c1 100644 --- a/geometry.py +++ b/geometry.py @@ -4,21 +4,23 @@ def rotation_matrix(alpha, rad=False): if not rad: - alpha /= 360/(2*np.pi) + alpha /= 360 / (2 * np.pi) return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]]) + def rotate(x, y, alpha, **kwargs): T = rotation_matrix(alpha, **kwargs) xy_rot = T @ np.vstack((x, y)) return xy_rot[0], xy_rot[1] + def plot_barriers(barriers, design, title=None, file_name=None): plt.figure() - alphas = np.linspace(0*np.pi,2*np.pi/4,100) + alphas = np.linspace(0 * np.pi, 2 * np.pi / 4, 100) for r in [design.rotor_r_min, design.rotor_r_max]: - plt.plot(r*np.cos(alphas), r*np.sin(alphas)) + plt.plot(r * np.cos(alphas), r * np.sin(alphas)) for barrier in barriers: - plt.plot(barrier[:,0], barrier[:,1]) + plt.plot(barrier[:, 0], barrier[:, 1]) plt.axis("equal") if title is not None: plt.title(title) @@ -26,15 +28,17 @@ def plot_barriers(barriers, design, title=None, file_name=None): plt.savefig(file_name) plt.close() + def analyze_results(Tor: np.ndarray) -> tuple[float, float, float]: TorAvg = np.mean(Tor[:-1]) - TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) - TorRippleRms = TorRmsAC/TorAvg*100 + TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1] - TorAvg))) + TorRippleRms = TorRmsAC / TorAvg * 100 return TorAvg, TorRmsAC, TorRippleRms - + + def print_results(TorAvg: float, TorRmsAC: float, TorRippleRms: float) -> None: print("\nTorque mean value: {:.2f} Nm".format(TorAvg)) # print("\nTorque ripple mean value: {:.2f} Nm".format(TorAvgAC)) print("\nTorque ripple rms value: {:.2f} Nm".format(TorRmsAC)) # print("\nTorque ripple relative value: {:.2f} %".format(TorRippleAvg)) - print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) \ No newline at end of file + print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) diff --git a/run.py b/run.py index 0c4ef63..0f0f24d 100644 --- a/run.py +++ b/run.py @@ -4,28 +4,36 @@ import numpy as np import pandas as pd from design import Design -from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta, save_params +from generators import ( + FourStupid, + HacklGenerator_OneLambda, + HacklGenerator_TwoLambdas, + HacklGenerator_OneLambdaTheta, + save_params, +) from geometry import analyze_results, plot_barriers project_name = "SynRM_test" design_name = "Design01" -path_data = os.path.join(os.getcwd(), 'data') -path_results = 'results' +path_data = os.path.join(os.getcwd(), "data") +path_results = "results" for path in [path_data, path_results]: os.makedirs(path, exist_ok=True) -file_name_aedt = f'{path_data}/{project_name}.aedt' +file_name_aedt = f"{path_data}/{project_name}.aedt" plot_design = True n_designs = 20 # Define constants AEDT_VERSION = "2024.1" NUM_CORES = 4 -NG_MODE = True #non-graphical mode -CLS_EXIT = True #close on exit +NG_MODE = True # non-graphical mode +CLS_EXIT = True # close on exit if not os.path.exists(file_name_aedt): design = Design.create( - project_name, design_name, file_name_aedt, + project_name, + design_name, + file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, new_desktop=False, @@ -65,7 +73,7 @@ Tor = design.compute(NUM_CORES) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan - else: + else: TorAvg, _, TorRippleRms = analyze_results(Tor) # Delete the rotor @@ -73,18 +81,20 @@ # Potentially save the design if plot_design: - title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %' - file_name = f'{path_results}/design_{generator.name}_{i}' - plot_barriers(barriers, design, title=title, file_name = f"{file_name}.png") + title = f"Torque mean value: {np.round(TorAvg, 2)} Nm, ripple relative value: {np.round(TorRippleRms, 2)} %" + file_name = f"{path_results}/design_{generator.name}_{i}" + plot_barriers(barriers, design, title=title, file_name=f"{file_name}.png") save_params(params, f"{file_name}.pkl") metadata_new = { - 'method': generator.name, - 'design': i, - 'T': TorAvg, - 'ripple': TorRippleRms, + "method": generator.name, + "design": i, + "T": TorAvg, + "ripple": TorRippleRms, } - metadata = pd.concat((metadata, pd.DataFrame([metadata_new])), ignore_index=True) + metadata = pd.concat( + (metadata, pd.DataFrame([metadata_new])), ignore_index=True + ) metadata.to_csv(f"{path_results}/metadata.csv", index=False) -design.close_project() \ No newline at end of file +design.close_project() diff --git a/run_optimization.py b/run_optimization.py index e98da44..9c53d06 100644 --- a/run_optimization.py +++ b/run_optimization.py @@ -6,9 +6,12 @@ import pandas as pd import pickle from design import Design -from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta -from geometry import analyze_results, plot_barriers +from generators import ( + FourStupid, +) +from geometry import analyze_results import torch + torch.set_default_dtype(torch.float64) from torch import Tensor from botorch.utils.transforms import normalize, unnormalize @@ -23,16 +26,19 @@ from botorch.utils.multi_objective.pareto import is_non_dominated from botorch.optim import optimize_acqf + def objective_transform(TorAvg, TorRippleRms): if pd.isnull(TorAvg): return -99999, -99999 else: - return TorAvg, -TorRippleRms/100 + return TorAvg, -TorRippleRms / 100 + def objective(Xs: Tensor, *args) -> Tensor: vals = [objective_single(X, *args) for X in Xs] return torch.stack(vals, dim=0) + def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: X = unnormalize(X, bounds) params = generator.X_to_params(X.numpy()) @@ -50,7 +56,7 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: Tor = design.compute(NUM_CORES) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan - else: + else: TorAvg, _, TorRippleRms = analyze_results(Tor) # Delete the rotor @@ -59,17 +65,21 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: f1, f2 = objective_transform(TorAvg, TorRippleRms) return torch.tensor([f1, f2]) + def init_points(root, method): results = pd.read_csv(f"{root}/metadata.csv") results = results[~results["T"].isnull()] results = results[results["method"] == method] - results["path"] = [f"results/design_{row['method']}_{row['design']}.pkl" for _, row in results.iterrows()] + results["path"] = [ + f"results/design_{row['method']}_{row['design']}.pkl" + for _, row in results.iterrows() + ] Xs, Ys = [], [] for _, r in results.iterrows(): - with open(r["path"], 'rb') as f: + with open(r["path"], "rb") as f: params = pickle.load(f) - + X = [] for x in params: if isinstance(x, Iterable): @@ -84,21 +94,24 @@ def init_points(root, method): return torch.stack(Xs), torch.stack(Ys) + project_name = "SynRM_test" design_name = "Design01" -path_data = os.path.join(os.getcwd(), 'data') +path_data = os.path.join(os.getcwd(), "data") os.makedirs(path_data, exist_ok=True) -file_name_aedt = f'{path_data}/{project_name}.aedt' +file_name_aedt = f"{path_data}/{project_name}.aedt" # Define constants AEDT_VERSION = "2024.1" NUM_CORES = 4 -NG_MODE = True #non-graphical mode -CLS_EXIT = True #close on exit +NG_MODE = True # non-graphical mode +CLS_EXIT = True # close on exit if not os.path.exists(file_name_aedt): design = Design.create( - project_name, design_name, file_name_aedt, + project_name, + design_name, + file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, new_desktop=False, @@ -126,11 +139,13 @@ def init_points(root, method): train_X, train_Y = init_points(root_init, method) train_X = normalize(train_X, bounds) bounds_normalized = normalize(bounds, bounds) -ref_point = torch.tensor([3.8,-0.3]) +ref_point = torch.tensor([3.8, -0.3]) + def objective_lambda(Xs): return objective(Xs, design, generator, bounds, NUM_CORES) + for _ in range(n_iters): # Fit surrogate model = SingleTaskGP(train_X, train_Y) @@ -139,10 +154,7 @@ def objective_lambda(Xs): # Compute Pareto front pareto_Y = train_Y[is_non_dominated(train_Y)] - partitioning = NondominatedPartitioning( - ref_point=ref_point, - Y=pareto_Y - ) + partitioning = NondominatedPartitioning(ref_point=ref_point, Y=pareto_Y) # Optimize acquisition function to select candidate points acq = qLogExpectedHypervolumeImprovement( @@ -162,6 +174,8 @@ def objective_lambda(Xs): new_Y = objective_lambda(candidates) train_X = torch.cat([train_X, candidates]) train_Y = torch.cat([train_Y, new_Y]) - np.savez(f"results_{method}.npz", train_X=unnormalize(train_X, bounds), train_Y=train_Y) + np.savez( + f"results_{method}.npz", train_X=unnormalize(train_X, bounds), train_Y=train_Y + ) -design.close_project() \ No newline at end of file +design.close_project() From 7aa66af5173122e4dd084bab5993806932425137 Mon Sep 17 00:00:00 2001 From: sadda Date: Mon, 9 Mar 2026 18:35:27 +0000 Subject: [PATCH 40/55] Revert "Ruff" This reverts commit 05e76ce1fcf878804bee5125045839fae766d284. --- design.py | 165 +++++++++++++++++--------------------------- generators.py | 153 ++++++++++++++++------------------------ geometry.py | 20 +++--- run.py | 44 +++++------- run_optimization.py | 52 +++++--------- 5 files changed, 167 insertions(+), 267 deletions(-) diff --git a/design.py b/design.py index 9cae6b8..819e974 100644 --- a/design.py +++ b/design.py @@ -1,22 +1,20 @@ +import os import numpy as np from ansys.aedt.core import Desktop, Maxwell2d from ansys.aedt.core.modeler.modeler_2d import Modeler2D - class Design: def __init__(self, m2d: Maxwell2d) -> None: self.set_parameters() self.m2d = m2d @classmethod - def create( - cls, project_name: str, design_name: str, file_name: str, **kwargs - ) -> "Design": + def create(cls, project_name: str, design_name: str, file_name: str, **kwargs) -> "Design": m2d = Maxwell2d( project=project_name, design=design_name, solution_type="TransientXY", - **kwargs, + **kwargs ) obj = cls(m2d) obj.create_stator() @@ -32,9 +30,9 @@ def load(cls, file_name: str, **kwargs) -> "Design": m2d.set_active_design("Design01") return cls(m2d) - + def set_parameters(self) -> None: - # materials + #materials self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" # main definitions @@ -45,7 +43,7 @@ def set_parameters(self) -> None: "SlotNumber": "36", "SlotType": "3", "DiaShaft": "25mm", - "StackLength": "85mm", + "StackLength": "85mm" } # stator slot self.slot_params = { @@ -65,57 +63,40 @@ def set_parameters(self) -> None: "CoilPitch": "9", # coil pitch in slots "SlotLiner": "0.3mm", "SpaceLayers": "0.2mm", - "Nc": "68", # turns per coil + "Nc": "68" # turns per coil } # model parameters PolePairs = 2 - f = 50 # [Hz] - RotSpeed = 60 * f / PolePairs # [rpm] + f = 50 #[Hz] + RotSpeed = 60*f/PolePairs #[rpm] self.mod_params = { "Poles": f"2*{PolePairs}", "ModelLength": "85mm", "SymmetryFactor": "Poles", "StatorSkewAngle": "0deg", } - # operation parameters + #operation parameters self.oper_params = { "Im": "1.5*sqrt(2)A", - "epsI": "pi/4", # current angle + "epsI": "pi/4", #current angle "InitPos": "-30deg", "f": f"{f}Hz", "RotSpeed": f"{RotSpeed}rpm", - "Nper": "1/6", # number of included periods - "PointPer": "101", # number of time points per period + "Nper": "1/6" + "", #number of included periods + "PointPer": "101" #number of time points per period } self.rot_points = [ - [ - "DiaShaft/2*cos(360deg/SymmetryFactor)", - "DiaShaft/2*sin(360deg/SymmetryFactor)", - "0mm", - ], - [ - "DiaShaft/2*cos(360deg/(2*SymmetryFactor))", - "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", - "0mm", - ], + ["DiaShaft/2*cos(360deg/SymmetryFactor)", "DiaShaft/2*sin(360deg/SymmetryFactor)", "0mm"], + ["DiaShaft/2*cos(360deg/(2*SymmetryFactor))", "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", "0mm"], ["DiaShaft/2", "0mm", "0mm"], ["DiaStatorGap/2-Airgap", "0mm", "0mm"], - [ - "(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", - "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", - "0mm", - ], - [ - "(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", - "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", - "0mm", - ], + ["(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", "0mm"], + ["(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", "0mm"], ] self.setup_name = "Setup1" - self.rotor_r_min = self.mm_to_str("geom_params", "DiaShaft") / 2 - self.rotor_r_max = self.mm_to_str( - "geom_params", "DiaStatorGap" - ) / 2 - self.mm_to_str("geom_params", "Airgap") + self.rotor_r_min = self.mm_to_str('geom_params', 'DiaShaft')/2 + self.rotor_r_max = self.mm_to_str('geom_params', 'DiaStatorGap')/2 - self.mm_to_str('geom_params', 'Airgap') def create_stator(self) -> None: m2d = self.m2d @@ -160,7 +141,7 @@ def create_stator(self) -> None: name="Shaft", ) - # motion setup + #motion setup m2d.assign_rotate_motion( assignment="Band", coordinate_system="Global", @@ -218,14 +199,10 @@ def create_stator(self) -> None: stator_id.color = (192, 192, 192) # rgb stator_id.transparency = 0.0 - # Winding + #Winding coil_id = modeler.create_rectangle( origin=["DiaStatorGap/2+Hs0+Hs1+SlotLiner", "-(Bs1/2-SlotLiner)", 0], - sizes=[ - "(Hs2+Rs/2-2*SlotLiner)/Layers-(Layers-1)*SpaceLayers/2", - "Bs1-2*SlotLiner", - 0, - ], + sizes=["(Hs2+Rs/2-2*SlotLiner)/Layers-(Layers-1)*SpaceLayers/2", "Bs1-2*SlotLiner", 0], name="Coil", material="copper", ) @@ -233,10 +210,7 @@ def create_stator(self) -> None: coil_id.transparency = 0.0 modeler.rotate(assignment=coil_id, axis="Z", angle="360deg/SlotNumber/2") coil_id.duplicate_around_axis( - axis="Z", - angle="360deg/SlotNumber", - clones="CoilPitch", - create_new_objects=True, + axis="Z", angle="360deg/SlotNumber", clones="CoilPitch", create_new_objects=True ) id_coils = modeler.get_objects_w_string(string_name="Coil", case_sensitive=True) @@ -286,14 +260,16 @@ def create_stator(self) -> None: ], assignment="Region", ) - m2d.assign_vector_potential(assignment=id_bc_az, vector_value=0, boundary="A0") + m2d.assign_vector_potential( + assignment=id_bc_az, vector_value=0, boundary="A0" + ) # Excitations I_A = "Im * cos(2*pi*f*time+epsI)" I_B = "Im * cos(2*pi*f*time-120deg+epsI)" I_C = "Im * cos(2*pi*f*time-240deg+epsI)" - # Define phase windings + #Define phase windings m2d.assign_coil( assignment=["Coil"], conductors_number="Nc", @@ -320,7 +296,9 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseA", ) - m2d.add_winding_coils(assignment="PhaseA", coils=["CS1", "CS2", "CS3"]) + m2d.add_winding_coils( + assignment="PhaseA", coils=["CS1", "CS2", "CS3"] + ) m2d.assign_coil( assignment=["Coil_6"], conductors_number="Nc", @@ -347,7 +325,9 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseB", ) - m2d.add_winding_coils(assignment="PhaseB", coils=["CS7", "CS8", "CS9"]) + m2d.add_winding_coils( + assignment="PhaseB", coils=["CS7", "CS8", "CS9"] + ) m2d.assign_coil( assignment=["Coil_3"], conductors_number="Nc", @@ -374,7 +354,9 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseC", ) - m2d.add_winding_coils(assignment="PhaseC", coils=["CS4", "CS5", "CS6"]) + m2d.add_winding_coils( + assignment="PhaseC", coils=["CS4", "CS5", "CS6"] + ) # Mesh operation m2d.mesh.assign_length_mesh( @@ -392,18 +374,18 @@ def create_stator(self) -> None: name="stator", ) - # core loss + #core loss m2d.set_core_losses("Stator", core_loss_on_field=False) - # inductance calculation + #inductance calculation m2d.change_inductance_computation( compute_transient_inductance=True, incremental_matrix=False ) - # model depth + #model depth m2d.model_depth = "StackLength" - # symmetry + #symmetry m2d.change_symmetry_multiplier("SymmetryFactor") - # Calculation setup + # Calculation setup setup = m2d.create_setup(name=self.setup_name) setup.props["StopTime"] = "Nper/f" setup.props["TimeStep"] = "1/(f*(PointPer-1))" @@ -414,7 +396,7 @@ def create_stator(self) -> None: setup.update() m2d.validate_simple() - # ooutput variables + #ooutput variables output_vars = { "pos": "(Moving1.Position -InitPos) * Poles/2", "cos0": "cos(pos)", @@ -442,34 +424,21 @@ def create_stator(self) -> None: for k, v in output_vars.items(): m2d.create_output_variable(k, v) - # Definitions for plots + #Definitions for plots post_params = { # reports - ( - "InducedVoltage(PhaseA)", - "InducedVoltage(PhaseB)", - "InducedVoltage(PhaseC)", - ): "InducedVoltage", + ("InducedVoltage(PhaseA)","InducedVoltage(PhaseB)","InducedVoltage(PhaseC)"): "InducedVoltage", ("Moving1.Torque"): "Torque", - ( - "InputCurrent(PhaseA)", - "InputCurrent(PhaseB)", - "InputCurrent(PhaseC)", - ): "Current", - ( - "FluxLinkage(PhaseA)", - "FluxLinkage(PhaseB)", - "FluxLinkage(PhaseC)", - ): "FluxLinkage", + ("InputCurrent(PhaseA)","InputCurrent(PhaseB)","InputCurrent(PhaseC)"): "Current", + ("FluxLinkage(PhaseA)","FluxLinkage(PhaseB)","FluxLinkage(PhaseC)",): "FluxLinkage", ("I_d", "I_q"): "Current_dq", ("Flux_d", "Flux_q"): "FluxLinkage_dq", ("Ui_d", "Ui_q"): "InducedVoltage_dq", ("L_d", "L_q"): "Inductance_dq", + } - # Create Report + #Create Report for k, v in post_params.items(): - expressions = ( - list(k) if isinstance(k, tuple) else [k] - ) # if multiple report, use list(k). Else, use k + expressions = list(k) if isinstance(k, tuple) else [k] #if multiple report, use list(k). Else, use k report = m2d.post.create_report( expressions=expressions, setup_sweep_name="", @@ -490,10 +459,7 @@ def add_rotor(self) -> None: assert isinstance(modeler, Modeler2D) rotor_id = modeler.create_polyline( - points=self.rot_points, - segment_type=["Arc", "Line", "Arc"], - cover_surface=True, - name="Rotor", + points=self.rot_points, segment_type=[ "Arc","Line", "Arc"], cover_surface=True, name="Rotor" ) self.rotor_id = rotor_id rotor_id.material_name = self.Fe @@ -507,19 +473,14 @@ def add_rotor_barrier(self, barrier_points, segment_type=None) -> None: # Round it for Ansys barrier_points = np.round(barrier_points, 6) - # Potentially add the z axis + # Potentially add the z axis if barrier_points.shape[1] == 2: - barrier_points = np.hstack( - (barrier_points, np.zeros((len(barrier_points), 1))) - ) - + barrier_points = np.hstack((barrier_points, np.zeros((len(barrier_points),1)))) + # Convert them into a string format and interpolate points_str = [[str(y) for y in x] for x in barrier_points] barrier_id = modeler.create_polyline( - points=points_str, - segment_type=segment_type, - cover_surface=True, - name="Barrier", + points=points_str, segment_type=segment_type, cover_surface=True, name="Barrier" ) # Remove the barrier @@ -538,7 +499,7 @@ def compute(self, NUM_CORES: int = 1): maximum_elements=None, name="rotor", ) - # core loss rotor + #core loss rotor m2d.set_core_losses("Rotor", core_loss_on_field=False) # Analyze @@ -551,14 +512,14 @@ def compute(self, NUM_CORES: int = 1): return solutions.data_magnitude() except AttributeError: return None - + def delete_rotor(self) -> None: assert isinstance(self.m2d.modeler, Modeler2D) self.m2d.modeler.delete(self.rotor_id) def save_design(self, file_name: str, **kwargs) -> None: - show = kwargs.pop("show", False) - view = kwargs.pop("view", "xy") + show = kwargs.pop('show', False) + view = kwargs.pop('view', 'xy') self.m2d.plot(show=show, output_file=file_name, view=view) def save_project(self, file_name: str | None = None) -> None: @@ -566,12 +527,12 @@ def save_project(self, file_name: str | None = None) -> None: self.m2d.save_project() else: self.m2d.save_project(file_name) - + def close_project(self) -> None: self.m2d.close_desktop() - + def mm_to_str(self, var, field) -> float: val = getattr(self, var)[field] - if not val.endswith("mm"): - raise Exception("val must end with mm") + if not val.endswith('mm'): + raise Exception('val must end with mm') return float(val[:-2]) diff --git a/generators.py b/generators.py index e1e5f32..1c1662f 100644 --- a/generators.py +++ b/generators.py @@ -15,24 +15,20 @@ def get_arc(R: float, start_deg: float, end_deg: float, n_points: int) -> np.nda return np.column_stack((x, y)) - def signed_distance(P: np.ndarray, a: float, b: float, c: float) -> np.ndarray: - return a * P[:, 0] + b * P[:, 1] + c - + return a * P[:,0] + b * P[:,1] + c def check_barrier(barrier: np.ndarray) -> None: if not np.array_equal(barrier[0], barrier[-1]): raise ValueError("Barrier must start and end at the same point") - def compute_r_max(design: Design, r_stator_end: float) -> float: dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") airgap = design.mm_to_str("geom_params", "Airgap") return (dia_stator_gap / 2.0) - airgap - r_stator_end - def save_params(params, file_name: str) -> None: - with open(file_name, "wb") as handle: + with open(file_name, 'wb') as handle: pickle.dump(params, handle, protocol=pickle.HIGHEST_PROTOCOL) @@ -41,7 +37,7 @@ def __init__(self, offset: float, n_curve: int = 500, n_flat: int = 20) -> None: self.offset = offset self.n_curve = n_curve self.n_flat = n_flat - + @abstractmethod def generate_barriers(self) -> list[np.ndarray]: pass @@ -68,7 +64,7 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: lam = np.linspace(0, 1, self.n_flat)[:, None] barriers = [] - for a, b, c in zip([1.0, -1.0], [-1.0, 1.0], [-self.offset, -self.offset]): + for a, b, c in zip([1., -1.], [-1., 1.], [-self.offset, -self.offset]): f = signed_distance(barrier, a, b, c) f_is_positive = f > 0 mask_positive = np.where(f_is_positive)[0] @@ -76,7 +72,7 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: barrier_new = [barrier[mask_positive[0]]] for i in range(len(mask_positive) - 1): i1 = mask_positive[i] - i2 = mask_positive[i + 1] + i2 = mask_positive[i+1] if i2 - i1 == 1: barrier_new.append(barrier[i2]) @@ -96,11 +92,8 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: return [x for barrier in barriers for x in self.split_barrier(barrier)] - class FourStupid(BarrierGenerator): - def __init__( - self, design, r_stator_end, der1=1.0, der2=1.0, symmetric=True, **kwargs - ) -> None: + def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwargs) -> None: self.R = compute_r_max(design, r_stator_end) self.design = design self.der1 = der1 @@ -110,31 +103,31 @@ def __init__( super().__init__(**kwargs) def _create_barrier( - self, - y_min, - w_min, - y_mid, - w_mid, - theta, - w_max, - ): - - theta1 = (theta + 45) / 180 * np.pi - x_max1 = self.R * np.cos(theta1) - y_max1 = self.R * np.sin(theta1) - - x1 = [0, x_max1 / 2, x_max1] + self, + y_min, + w_min, + y_mid, + w_mid, + theta, + w_max, + ): + + theta1 = (theta+45) / 180*np.pi + x_max1 = self.R*np.cos(theta1) + y_max1 = self.R*np.sin(theta1) + + x1 = [0, x_max1/2, x_max1] y1 = [y_min, y_mid, y_max1] - f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1, self.der1))) + f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,self.der1))) theta2 = theta1 + w_max / self.design.rotor_r_max - x_max2 = self.R * np.cos(theta2) - y_max2 = self.R * np.sin(theta2) + x_max2 = self.R*np.cos(theta2) + y_max2 = self.R*np.sin(theta2) - s = w_mid / np.sqrt(1 + self.der1**2 / 4) - x2 = [0, x1[1] - s * self.der1 / 2, x_max2] - y2 = [y_min + w_min, y1[1] + s, y_max2] - f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1, self.der2))) + s = w_mid / np.sqrt(1 + self.der1**2/4) + x2 = [0, x1[1]-s*self.der1/2, x_max2] + y2 = [y_min+w_min, y1[1]+s, y_max2] + f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,self.der2))) x_interp1 = np.linspace(x1[0], x1[-1], self.n_curve) x_interp2 = np.linspace(x2[0], x2[-1], self.n_curve) @@ -154,20 +147,20 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: return lb, ub def random_parameters(self): - rand1 = 0.5 * np.random.random(4) - rand2 = 13 + 4 * np.random.random() - rand3 = 2 + 3 * np.random.random() - rand4 = 2 + 3 * np.random.random() - rand5 = 2 + 2 * np.random.random() + rand1 = 0.5*np.random.random(4) + rand2 = 13+4*np.random.random() + rand3 = 2+3*np.random.random() + rand4 = 2+3*np.random.random() + rand5 = 2+2*np.random.random() return rand1, rand2, rand3, rand4, rand5 def set_parameters(self, params) -> None: rand1, rand2, rand3, rand4, rand5 = params self.w_mins = self.w_mins_base + rand1 y_min0 = rand2 - y_min1 = y_min0 + self.w_mins[0] + rand3 - y_min2 = y_min1 + self.w_mins[1] + rand4 - y_min3 = y_min2 + self.w_mins[2] + rand5 + y_min1 = y_min0+self.w_mins[0] + rand3 + y_min2 = y_min1+self.w_mins[1] + rand4 + y_min3 = y_min2+self.w_mins[2] + rand5 self.w_mids = self.w_mins self.y_mins = np.array([y_min0, y_min1, y_min2, y_min3]) self.y_mids = self.y_mins + np.array([2, 1.5, 1, 0.5]) @@ -179,38 +172,30 @@ def X_to_params(self, X: np.ndarray): def generate_barriers(self) -> list[np.ndarray]: barriers = [] - for args in zip( - self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs - ): + for args in zip(self.y_mins, self.w_mins, self.y_mids, self.w_mids, self.thetas, self.w_maxs): x_all, y_all = self._create_barrier(*args) xy_all = np.vstack((x_all, y_all)).T barriers.append(xy_all) return barriers - + class AbstractHacklGenerator(BarrierGenerator): def __init__(self, design, r_stator_end, **kwargs): self.R = compute_r_max(design, r_stator_end) - self.phis_inner_min = np.asarray([4.0, 17.3, 30.6]) + self.phis_inner_min = np.asarray([4., 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) - self.phis_outer_min = np.asarray([10.6, 24.0, 37.3]) - self.phis_outer_max = np.asarray([16.8, 30.2, 44.0]) + self.phis_outer_min = np.asarray([10.6, 24., 37.3]) + self.phis_outer_max = np.asarray([16.8, 30.2, 44.]) self.n_barriers = len(self.phis_inner_min) super().__init__(**kwargs) @abstractmethod - def _bezier_point( - self, x: float, y: float, is_inner: bool, order: int - ) -> tuple[float, float]: + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: pass def random_parameters(self): - phis_inner = self.phis_inner_min + ( - self.phis_inner_max - self.phis_inner_min - ) * np.random.rand(self.n_barriers) - phis_outer = self.phis_outer_min + ( - self.phis_outer_max - self.phis_outer_min - ) * np.random.rand(self.n_barriers) + phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) + phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) return phis_inner, phis_outer def set_parameters(self, params) -> None: @@ -220,21 +205,17 @@ def set_parameters(self, params) -> None: def generate_barriers(self) -> list[np.ndarray]: barriers = [] - for i, (phi_inner, phi_outer) in enumerate( - zip(self.phis_inner, self.phis_outer) - ): + for i, (phi_inner, phi_outer) in enumerate(zip(self.phis_inner, self.phis_outer)): # Generate the longer part of the barrier pts_outer = self._get_bezier_curve(phi_outer, False, i) pts_inner = self._get_bezier_curve(phi_inner, True, i) - + # Generate the arcs at the end arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner, self.n_flat) arc_bottom = get_arc(self.R, phi_inner, phi_outer, self.n_flat) - + # Merge them together - barrier = np.concatenate( - (pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:]) - ) + barrier = np.concatenate((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) barriers.append(barrier) return barriers @@ -252,8 +233,8 @@ def _get_bezier_curve(self, phi_deg, is_inner, i): # Formula 33: Polynomial expansion to generate continuous points z_vals = np.linspace(0, 1, self.n_curve)[:, None] return ( - (1 - z_vals) ** 3 * r0 - + 3 * (1 - z_vals) ** 2 * z_vals * r1 + (1 - z_vals)**3 * r0 + + 3 * (1 - z_vals)**2 * z_vals * r1 + 3 * (1 - z_vals) * z_vals**2 * r2 + z_vals**3 * r3 ) @@ -270,12 +251,12 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: lb = np.concatenate((self.phis_inner_min, self.phis_outer_min, [self.lam_min])) ub = np.concatenate((self.phis_inner_max, self.phis_outer_max, [self.lam_max])) return lb, ub - + def random_parameters(self): pars = super().random_parameters() - lam = self.lam_min + (self.lam_max - self.lam_min) * np.random.rand(1)[0] + lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] return *pars, lam - + def set_parameters(self, params) -> None: phis_inner, phis_outer, lam = params super().set_parameters((phis_inner, phis_outer)) @@ -283,10 +264,8 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6] - - def _bezier_point( - self, x: float, y: float, is_inner: bool, order: int - ) -> tuple[float, float]: + + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y @@ -298,9 +277,7 @@ def __init__(self, design, r_stator_end, **kwargs): def random_parameters(self): pars = super().random_parameters() - theta = ( - self.theta_min + (self.theta_max - self.theta_min) * np.random.rand(1)[0] - ) + theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] return *pars, theta def set_parameters(self, params) -> None: @@ -310,10 +287,8 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6], X[7] - - def _bezier_point( - self, x: float, y: float, is_inner: bool, order: int - ) -> tuple[float, float]: + + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: x, y = super()._bezier_point(x, y, is_inner, order) if order < 2: xx, yy = rotate(x, y, self.theta) @@ -331,14 +306,8 @@ def __init__(self, design, r_stator_end, **kwargs): def random_parameters(self): pars = super().random_parameters() - lam_inner = ( - self.lam_inner_min - + (self.lam_inner_max - self.lam_inner_min) * np.random.rand(1)[0] - ) - lam_outer = ( - self.lam_outer_min - + (self.lam_outer_max - self.lam_outer_min) * np.random.rand(1)[0] - ) + lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] + lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] return *pars, lam_inner, lam_outer def set_parameters(self, params) -> None: @@ -350,8 +319,6 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6], X[7] - def _bezier_point( - self, x: float, y: float, is_inner: bool, order: int - ) -> tuple[float, float]: + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: lam = self.lam_inner if is_inner else self.lam_outer return lam * x + (1 - lam) * y, y diff --git a/geometry.py b/geometry.py index 3e659c1..328ba07 100644 --- a/geometry.py +++ b/geometry.py @@ -4,23 +4,21 @@ def rotation_matrix(alpha, rad=False): if not rad: - alpha /= 360 / (2 * np.pi) + alpha /= 360/(2*np.pi) return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]]) - def rotate(x, y, alpha, **kwargs): T = rotation_matrix(alpha, **kwargs) xy_rot = T @ np.vstack((x, y)) return xy_rot[0], xy_rot[1] - def plot_barriers(barriers, design, title=None, file_name=None): plt.figure() - alphas = np.linspace(0 * np.pi, 2 * np.pi / 4, 100) + alphas = np.linspace(0*np.pi,2*np.pi/4,100) for r in [design.rotor_r_min, design.rotor_r_max]: - plt.plot(r * np.cos(alphas), r * np.sin(alphas)) + plt.plot(r*np.cos(alphas), r*np.sin(alphas)) for barrier in barriers: - plt.plot(barrier[:, 0], barrier[:, 1]) + plt.plot(barrier[:,0], barrier[:,1]) plt.axis("equal") if title is not None: plt.title(title) @@ -28,17 +26,15 @@ def plot_barriers(barriers, design, title=None, file_name=None): plt.savefig(file_name) plt.close() - def analyze_results(Tor: np.ndarray) -> tuple[float, float, float]: TorAvg = np.mean(Tor[:-1]) - TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1] - TorAvg))) - TorRippleRms = TorRmsAC / TorAvg * 100 + TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) + TorRippleRms = TorRmsAC/TorAvg*100 return TorAvg, TorRmsAC, TorRippleRms - - + def print_results(TorAvg: float, TorRmsAC: float, TorRippleRms: float) -> None: print("\nTorque mean value: {:.2f} Nm".format(TorAvg)) # print("\nTorque ripple mean value: {:.2f} Nm".format(TorAvgAC)) print("\nTorque ripple rms value: {:.2f} Nm".format(TorRmsAC)) # print("\nTorque ripple relative value: {:.2f} %".format(TorRippleAvg)) - print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) + print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) \ No newline at end of file diff --git a/run.py b/run.py index 0f0f24d..0c4ef63 100644 --- a/run.py +++ b/run.py @@ -4,36 +4,28 @@ import numpy as np import pandas as pd from design import Design -from generators import ( - FourStupid, - HacklGenerator_OneLambda, - HacklGenerator_TwoLambdas, - HacklGenerator_OneLambdaTheta, - save_params, -) +from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta, save_params from geometry import analyze_results, plot_barriers project_name = "SynRM_test" design_name = "Design01" -path_data = os.path.join(os.getcwd(), "data") -path_results = "results" +path_data = os.path.join(os.getcwd(), 'data') +path_results = 'results' for path in [path_data, path_results]: os.makedirs(path, exist_ok=True) -file_name_aedt = f"{path_data}/{project_name}.aedt" +file_name_aedt = f'{path_data}/{project_name}.aedt' plot_design = True n_designs = 20 # Define constants AEDT_VERSION = "2024.1" NUM_CORES = 4 -NG_MODE = True # non-graphical mode -CLS_EXIT = True # close on exit +NG_MODE = True #non-graphical mode +CLS_EXIT = True #close on exit if not os.path.exists(file_name_aedt): design = Design.create( - project_name, - design_name, - file_name_aedt, + project_name, design_name, file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, new_desktop=False, @@ -73,7 +65,7 @@ Tor = design.compute(NUM_CORES) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan - else: + else: TorAvg, _, TorRippleRms = analyze_results(Tor) # Delete the rotor @@ -81,20 +73,18 @@ # Potentially save the design if plot_design: - title = f"Torque mean value: {np.round(TorAvg, 2)} Nm, ripple relative value: {np.round(TorRippleRms, 2)} %" - file_name = f"{path_results}/design_{generator.name}_{i}" - plot_barriers(barriers, design, title=title, file_name=f"{file_name}.png") + title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %' + file_name = f'{path_results}/design_{generator.name}_{i}' + plot_barriers(barriers, design, title=title, file_name = f"{file_name}.png") save_params(params, f"{file_name}.pkl") metadata_new = { - "method": generator.name, - "design": i, - "T": TorAvg, - "ripple": TorRippleRms, + 'method': generator.name, + 'design': i, + 'T': TorAvg, + 'ripple': TorRippleRms, } - metadata = pd.concat( - (metadata, pd.DataFrame([metadata_new])), ignore_index=True - ) + metadata = pd.concat((metadata, pd.DataFrame([metadata_new])), ignore_index=True) metadata.to_csv(f"{path_results}/metadata.csv", index=False) -design.close_project() +design.close_project() \ No newline at end of file diff --git a/run_optimization.py b/run_optimization.py index 9c53d06..e98da44 100644 --- a/run_optimization.py +++ b/run_optimization.py @@ -6,12 +6,9 @@ import pandas as pd import pickle from design import Design -from generators import ( - FourStupid, -) -from geometry import analyze_results +from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta +from geometry import analyze_results, plot_barriers import torch - torch.set_default_dtype(torch.float64) from torch import Tensor from botorch.utils.transforms import normalize, unnormalize @@ -26,19 +23,16 @@ from botorch.utils.multi_objective.pareto import is_non_dominated from botorch.optim import optimize_acqf - def objective_transform(TorAvg, TorRippleRms): if pd.isnull(TorAvg): return -99999, -99999 else: - return TorAvg, -TorRippleRms / 100 - + return TorAvg, -TorRippleRms/100 def objective(Xs: Tensor, *args) -> Tensor: vals = [objective_single(X, *args) for X in Xs] return torch.stack(vals, dim=0) - def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: X = unnormalize(X, bounds) params = generator.X_to_params(X.numpy()) @@ -56,7 +50,7 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: Tor = design.compute(NUM_CORES) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan - else: + else: TorAvg, _, TorRippleRms = analyze_results(Tor) # Delete the rotor @@ -65,21 +59,17 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: f1, f2 = objective_transform(TorAvg, TorRippleRms) return torch.tensor([f1, f2]) - def init_points(root, method): results = pd.read_csv(f"{root}/metadata.csv") results = results[~results["T"].isnull()] results = results[results["method"] == method] - results["path"] = [ - f"results/design_{row['method']}_{row['design']}.pkl" - for _, row in results.iterrows() - ] + results["path"] = [f"results/design_{row['method']}_{row['design']}.pkl" for _, row in results.iterrows()] Xs, Ys = [], [] for _, r in results.iterrows(): - with open(r["path"], "rb") as f: + with open(r["path"], 'rb') as f: params = pickle.load(f) - + X = [] for x in params: if isinstance(x, Iterable): @@ -94,24 +84,21 @@ def init_points(root, method): return torch.stack(Xs), torch.stack(Ys) - project_name = "SynRM_test" design_name = "Design01" -path_data = os.path.join(os.getcwd(), "data") +path_data = os.path.join(os.getcwd(), 'data') os.makedirs(path_data, exist_ok=True) -file_name_aedt = f"{path_data}/{project_name}.aedt" +file_name_aedt = f'{path_data}/{project_name}.aedt' # Define constants AEDT_VERSION = "2024.1" NUM_CORES = 4 -NG_MODE = True # non-graphical mode -CLS_EXIT = True # close on exit +NG_MODE = True #non-graphical mode +CLS_EXIT = True #close on exit if not os.path.exists(file_name_aedt): design = Design.create( - project_name, - design_name, - file_name_aedt, + project_name, design_name, file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, new_desktop=False, @@ -139,13 +126,11 @@ def init_points(root, method): train_X, train_Y = init_points(root_init, method) train_X = normalize(train_X, bounds) bounds_normalized = normalize(bounds, bounds) -ref_point = torch.tensor([3.8, -0.3]) - +ref_point = torch.tensor([3.8,-0.3]) def objective_lambda(Xs): return objective(Xs, design, generator, bounds, NUM_CORES) - for _ in range(n_iters): # Fit surrogate model = SingleTaskGP(train_X, train_Y) @@ -154,7 +139,10 @@ def objective_lambda(Xs): # Compute Pareto front pareto_Y = train_Y[is_non_dominated(train_Y)] - partitioning = NondominatedPartitioning(ref_point=ref_point, Y=pareto_Y) + partitioning = NondominatedPartitioning( + ref_point=ref_point, + Y=pareto_Y + ) # Optimize acquisition function to select candidate points acq = qLogExpectedHypervolumeImprovement( @@ -174,8 +162,6 @@ def objective_lambda(Xs): new_Y = objective_lambda(candidates) train_X = torch.cat([train_X, candidates]) train_Y = torch.cat([train_Y, new_Y]) - np.savez( - f"results_{method}.npz", train_X=unnormalize(train_X, bounds), train_Y=train_Y - ) + np.savez(f"results_{method}.npz", train_X=unnormalize(train_X, bounds), train_Y=train_Y) -design.close_project() +design.close_project() \ No newline at end of file From 7401dffa6223a4b3c2f3e7bd0b860a177c1de8c5 Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 10 Mar 2026 08:47:58 +0000 Subject: [PATCH 41/55] Added bounds --- generators.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/generators.py b/generators.py index 1c1662f..cef310e 100644 --- a/generators.py +++ b/generators.py @@ -275,6 +275,12 @@ def __init__(self, design, r_stator_end, **kwargs): self.theta_max = 20 super().__init__(design, r_stator_end, **kwargs) + @property + def bounds(self) -> tuple[np.ndarray, np.ndarray]: + lb = np.concatenate((self.phis_inner_min, self.phis_outer_min, [self.lam_min, self.theta_min])) + ub = np.concatenate((self.phis_inner_max, self.phis_outer_max, [self.lam_max, self.theta_max])) + return lb, ub + def random_parameters(self): pars = super().random_parameters() theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] @@ -304,6 +310,12 @@ def __init__(self, design, r_stator_end, **kwargs): self.lam_outer_max = 0.45 super().__init__(design, r_stator_end, **kwargs) + @property + def bounds(self) -> tuple[np.ndarray, np.ndarray]: + lb = np.concatenate((self.phis_inner_min, self.phis_outer_min, [self.lam_inner_min, self.lam_outer_min])) + ub = np.concatenate((self.phis_inner_max, self.phis_outer_max, [self.lam_inner_max, self.lam_outer_max])) + return lb, ub + def random_parameters(self): pars = super().random_parameters() lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] From 3c084b287e5bc73a494a1fff68bc0a757a8d5adc Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 10 Mar 2026 09:49:01 +0000 Subject: [PATCH 42/55] Added barriers check --- generators.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/generators.py b/generators.py index cef310e..4ac93b5 100644 --- a/generators.py +++ b/generators.py @@ -2,6 +2,7 @@ import numpy as np from scipy.interpolate import CubicSpline import pickle +from shapely.geometry import LineString from design import Design from geometry import rotate from typing import Any @@ -23,9 +24,10 @@ def check_barrier(barrier: np.ndarray) -> None: raise ValueError("Barrier must start and end at the same point") def compute_r_max(design: Design, r_stator_end: float) -> float: - dia_stator_gap = design.mm_to_str("geom_params", "DiaStatorGap") - airgap = design.mm_to_str("geom_params", "Airgap") - return (dia_stator_gap / 2.0) - airgap - r_stator_end + return design.rotor_r_max - r_stator_end + +def compute_r_min(design: Design) -> float: + return design.rotor_r_min def save_params(params, file_name: str) -> None: with open(file_name, 'wb') as handle: @@ -58,6 +60,22 @@ def X_to_params(self, X: np.ndarray) -> tuple[Any, ...]: def name(self) -> str: return self.__class__.__name__ + def feasible_barriers(self, barriers: list[np.ndarray]) -> bool: + n = len(barriers) + curves = [LineString(barrier) for barrier in barriers] + + # Check if curve intersects itself + for curve in curves: + if not curve.is_simple: + return False + + # Check if two curves intersect themselves + for i in range(n): + for j in range(i+1,n): + if curves[i].intersects(curves[j]): + return False + return True + def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: assert self.offset is not None check_barrier(barrier) From 913c8d6bad9cf25b25f172a6208aae12137fc74e Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 10 Mar 2026 09:55:30 +0000 Subject: [PATCH 43/55] Better saving diameters --- generators.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/generators.py b/generators.py index 4ac93b5..6e666f6 100644 --- a/generators.py +++ b/generators.py @@ -23,19 +23,17 @@ def check_barrier(barrier: np.ndarray) -> None: if not np.array_equal(barrier[0], barrier[-1]): raise ValueError("Barrier must start and end at the same point") -def compute_r_max(design: Design, r_stator_end: float) -> float: - return design.rotor_r_max - r_stator_end - -def compute_r_min(design: Design) -> float: - return design.rotor_r_min - def save_params(params, file_name: str) -> None: with open(file_name, 'wb') as handle: pickle.dump(params, handle, protocol=pickle.HIGHEST_PROTOCOL) class BarrierGenerator(ABC): - def __init__(self, offset: float, n_curve: int = 500, n_flat: int = 20) -> None: + # TODO: rename offset + def __init__(self, design: Design, r_stator_end: float, offset: None | float = None, n_curve: int = 500, n_flat: int = 20) -> None: + self.r_max = design.rotor_r_max + self.r_min = design.rotor_r_min + self.R = self.r_max - r_stator_end self.offset = offset self.n_curve = n_curve self.n_flat = n_flat @@ -112,13 +110,11 @@ def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: class FourStupid(BarrierGenerator): def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwargs) -> None: - self.R = compute_r_max(design, r_stator_end) - self.design = design self.der1 = der1 self.der2 = der2 self.w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0 self.symmetric = symmetric - super().__init__(**kwargs) + super().__init__(design, r_stator_end, **kwargs) def _create_barrier( self, @@ -138,7 +134,7 @@ def _create_barrier( y1 = [y_min, y_mid, y_max1] f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,self.der1))) - theta2 = theta1 + w_max / self.design.rotor_r_max + theta2 = theta1 + w_max / self.r_max x_max2 = self.R*np.cos(theta2) y_max2 = self.R*np.sin(theta2) @@ -199,13 +195,12 @@ def generate_barriers(self) -> list[np.ndarray]: class AbstractHacklGenerator(BarrierGenerator): def __init__(self, design, r_stator_end, **kwargs): - self.R = compute_r_max(design, r_stator_end) self.phis_inner_min = np.asarray([4., 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) self.phis_outer_min = np.asarray([10.6, 24., 37.3]) self.phis_outer_max = np.asarray([16.8, 30.2, 44.]) self.n_barriers = len(self.phis_inner_min) - super().__init__(**kwargs) + super().__init__(design, r_stator_end, **kwargs) @abstractmethod def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: From e15bfee96a21c12621f73d0d9d7f7fa0699345dc Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 10 Mar 2026 10:06:08 +0000 Subject: [PATCH 44/55] Added checks on diameters --- generators.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/generators.py b/generators.py index 6e666f6..c4ae637 100644 --- a/generators.py +++ b/generators.py @@ -58,7 +58,13 @@ def X_to_params(self, X: np.ndarray) -> tuple[Any, ...]: def name(self) -> str: return self.__class__.__name__ - def feasible_barriers(self, barriers: list[np.ndarray]) -> bool: + def feasible_barriers(self, barriers: list[np.ndarray], tol: float = 1e-8) -> bool: + # Check whether the barriers are within bounds + for barrier in barriers: + norm_points = np.linalg.norm(barrier, axis=1) + if any(norm_points > self.R + tol) or any(norm_points < self.r_min - tol): + return False + n = len(barriers) curves = [LineString(barrier) for barrier in barriers] From cb4066d044837d7de44beb8e3bc2b2e3fcd1cd5d Mon Sep 17 00:00:00 2001 From: sadda Date: Wed, 11 Mar 2026 04:34:16 +0000 Subject: [PATCH 45/55] Clean up --- generate_designs.ipynb | 117 ----------------------------------------- 1 file changed, 117 deletions(-) delete mode 100644 generate_designs.ipynb diff --git a/generate_designs.ipynb b/generate_designs.ipynb deleted file mode 100644 index 7db8ebd..0000000 --- a/generate_designs.ipynb +++ /dev/null @@ -1,117 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "import numpy as np\n", - "from design import Design\n", - "from generators import FourStupid, HacklGenerator\n", - "from geometry import analyze_results, plot_barriers\n", - "\n", - "project_name = \"SynRM_test\"\n", - "design_name = \"Design01\"\n", - "path_data = os.path.join(os.getcwd(), 'data')\n", - "path_results = 'results'\n", - "for path in [path_data, path_results]:\n", - " os.makedirs(path, exist_ok=True)\n", - "file_name_aedt = f'{path_data}/{project_name}.aedt'\n", - "plot_design = True\n", - "n_designs = 10\n", - "\n", - "# Define constants\n", - "AEDT_VERSION = \"2024.1\"\n", - "NUM_CORES = 4\n", - "NG_MODE = True #non-graphical mode\n", - "CLS_EXIT = True #close on exit\n", - "\n", - "if not os.path.exists(file_name_aedt):\n", - " design = Design.create(\n", - " project_name, design_name, file_name_aedt,\n", - " version=AEDT_VERSION,\n", - " non_graphical=NG_MODE,\n", - " new_desktop=False,\n", - " close_on_exit=CLS_EXIT,\n", - " )\n", - "else:\n", - " design = Design.load(\n", - " file_name_aedt,\n", - " version=AEDT_VERSION,\n", - " non_graphical=NG_MODE,\n", - " new_desktop=False,\n", - " close_on_exit=CLS_EXIT,\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "r_stator_end = 0.7\n", - "offset = 0.7 / 2\n", - "generator_stupid = FourStupid(design, r_stator_end, offset=offset)\n", - "generator_hackl = HacklGenerator(design, r_stator_end, offset=offset)\n", - "generators = [generator_stupid, generator_hackl]\n", - "\n", - "for i in range(0, n_designs):\n", - " for generator in generators:\n", - " barriers1 = generator.random_barriers()\n", - " barriers2 = generator.split_barriers(barriers1)\n", - " for j, barriers in enumerate([barriers1, barriers2]):\n", - " design.add_rotor()\n", - "\n", - " for barrier in barriers:\n", - " design.add_rotor_barrier(barrier)\n", - "\n", - " # Compute the torque\n", - " Tor = design.compute(NUM_CORES)\n", - " TorAvg, _, TorRippleRms = analyze_results(Tor)\n", - "\n", - " # Delete the rotor\n", - " design.delete_rotor()\n", - "\n", - " # Potentially save the design\n", - " if plot_design:\n", - " title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %'\n", - " file_name = f'{path_results}/design_{generator.name}_{i}_{j}'\n", - " plot_barriers(barriers, design, title=title, file_name = f\"{file_name}.png\")\n", - " generator.save_barriers(f\"{file_name}.npz\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "design.close_project()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "venv", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.11" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -} From 458b2fefee0b2893e5e28a3eb87f47c41d214146 Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:45:13 +0000 Subject: [PATCH 46/55] Bounds changed --- generators.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/generators.py b/generators.py index c4ae637..f24cb32 100644 --- a/generators.py +++ b/generators.py @@ -323,10 +323,10 @@ def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple class HacklGenerator_TwoLambdas(AbstractHacklGenerator): def __init__(self, design, r_stator_end, **kwargs): - self.lam_inner_min = 0.25 - self.lam_inner_max = 0.45 - self.lam_outer_min = 0.25 - self.lam_outer_max = 0.45 + self.lam_inner_min = 0.2 + self.lam_inner_max = 0.5 + self.lam_outer_min = 0.2 + self.lam_outer_max = 0.5 super().__init__(design, r_stator_end, **kwargs) @property From e343c9e782b9d0a7b60043e824fe5505e7b377d3 Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:46:18 +0000 Subject: [PATCH 47/55] Clean up --- geometry.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/geometry.py b/geometry.py index 328ba07..b67ec07 100644 --- a/geometry.py +++ b/geometry.py @@ -31,10 +31,3 @@ def analyze_results(Tor: np.ndarray) -> tuple[float, float, float]: TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) TorRippleRms = TorRmsAC/TorAvg*100 return TorAvg, TorRmsAC, TorRippleRms - -def print_results(TorAvg: float, TorRmsAC: float, TorRippleRms: float) -> None: - print("\nTorque mean value: {:.2f} Nm".format(TorAvg)) - # print("\nTorque ripple mean value: {:.2f} Nm".format(TorAvgAC)) - print("\nTorque ripple rms value: {:.2f} Nm".format(TorRmsAC)) - # print("\nTorque ripple relative value: {:.2f} %".format(TorRippleAvg)) - print("\nTorque ripple relative value: {:.2f} %\n".format(TorRippleRms)) \ No newline at end of file From eea815a614aa20018952090d1ba27e5f48189ccd Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:48:49 +0000 Subject: [PATCH 48/55] Refactor --- machine_design/__init__.py | 3 +++ design.py => machine_design/design.py | 0 generators.py => machine_design/generators.py | 0 geometry.py => machine_design/geometry.py | 0 4 files changed, 3 insertions(+) create mode 100644 machine_design/__init__.py rename design.py => machine_design/design.py (100%) rename generators.py => machine_design/generators.py (100%) rename geometry.py => machine_design/geometry.py (100%) diff --git a/machine_design/__init__.py b/machine_design/__init__.py new file mode 100644 index 0000000..3839c2c --- /dev/null +++ b/machine_design/__init__.py @@ -0,0 +1,3 @@ +from .design import Design +from .generators import HacklGenerator_OneLambda, HacklGenerator_OneLambdaTheta, HacklGenerator_TwoLambdas, FourStupid, save_params +from .geometry import analyze_results, plot_barriers, rotate \ No newline at end of file diff --git a/design.py b/machine_design/design.py similarity index 100% rename from design.py rename to machine_design/design.py diff --git a/generators.py b/machine_design/generators.py similarity index 100% rename from generators.py rename to machine_design/generators.py diff --git a/geometry.py b/machine_design/geometry.py similarity index 100% rename from geometry.py rename to machine_design/geometry.py From ac438adce0bb92c88ee9e2051feb9fd71fb69c2b Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:49:44 +0000 Subject: [PATCH 49/55] Added constraints --- run_optimization.py | 48 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 37 insertions(+), 11 deletions(-) diff --git a/run_optimization.py b/run_optimization.py index e98da44..ca266af 100644 --- a/run_optimization.py +++ b/run_optimization.py @@ -113,12 +113,14 @@ def init_points(root, method): close_on_exit=CLS_EXIT, ) +max_ripple = 0.1 n_iters = 40 batch_size = 4 r_stator_end = 0.7 offset = 0.7 / 2 -# generator = HacklGenerator_OneLambda(design, r_stator_end, offset=offset) -generator = FourStupid(design, r_stator_end, offset=offset) +generator = HacklGenerator_OneLambda(design, r_stator_end, offset=offset) +# generator = FourStupid(design, r_stator_end, offset=offset) +# generator = HacklGenerator_TwoLambdas(design, r_stator_end, offset=offset) bounds = torch.from_numpy(np.vstack(generator.bounds)) root_init = "results" @@ -126,11 +128,15 @@ def init_points(root, method): train_X, train_Y = init_points(root_init, method) train_X = normalize(train_X, bounds) bounds_normalized = normalize(bounds, bounds) -ref_point = torch.tensor([3.8,-0.3]) +ref_point = torch.tensor([3.8,-max_ripple]) def objective_lambda(Xs): return objective(Xs, design, generator, bounds, NUM_CORES) +def ripple_constraint(Y): + ripple = -Y[...,1] + return ripple - max_ripple + for _ in range(n_iters): # Fit surrogate model = SingleTaskGP(train_X, train_Y) @@ -144,19 +150,39 @@ def objective_lambda(Xs): Y=pareto_Y ) - # Optimize acquisition function to select candidate points + # Define acquisition function acq = qLogExpectedHypervolumeImprovement( model=model, ref_point=ref_point.tolist(), partitioning=partitioning, + constraints=[ripple_constraint], ) - candidates, _ = optimize_acqf( - acq_function=acq, - bounds=bounds_normalized, - q=batch_size, - num_restarts=10, - raw_samples=128, - ) + + # Optimize acquisition function to select candidate points. Reject unfeasible points + candidates_feasible = [] + while True: + candidates, _ = optimize_acqf( + acq_function=acq, + bounds=bounds_normalized, + q=batch_size, + num_restarts=10, + raw_samples=128, + ) + for candidate in candidates: + candidate_normalized = unnormalize(candidate, bounds) + params = generator.X_to_params(candidate_normalized.numpy()) + + generator.set_parameters(params) + barriers = generator.generate_barriers() + barriers = generator.split_barriers(barriers) + feasible = generator.feasible_barriers(barriers) + if feasible: + candidates_feasible.append(candidate) + if len(candidates_feasible) >= batch_size: + break + if len(candidates_feasible) >= batch_size: + break + candidates = torch.stack(candidates_feasible) # Evaluate candidates new_Y = objective_lambda(candidates) From 3053a8aea57a88912bd61c3f97f7e89d17766c3a Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:50:25 +0000 Subject: [PATCH 50/55] Refactor --- run.py | 5 ++--- run_optimization.py | 5 ++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/run.py b/run.py index 0c4ef63..1ca3622 100644 --- a/run.py +++ b/run.py @@ -3,9 +3,8 @@ import os import numpy as np import pandas as pd -from design import Design -from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta, save_params -from geometry import analyze_results, plot_barriers +from machine_design import Design, analyze_results, plot_barriers +from machine_design import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta, save_params project_name = "SynRM_test" design_name = "Design01" diff --git a/run_optimization.py b/run_optimization.py index ca266af..21b22b3 100644 --- a/run_optimization.py +++ b/run_optimization.py @@ -5,9 +5,8 @@ import numpy as np import pandas as pd import pickle -from design import Design -from generators import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta -from geometry import analyze_results, plot_barriers +from machine_design import Design, analyze_results, plot_barriers +from machine_design import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta import torch torch.set_default_dtype(torch.float64) from torch import Tensor From 8254b1659497a4754961509f724c9a50c057314a Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:53:14 +0000 Subject: [PATCH 51/55] Added pyproject --- pyproject.toml | 51 ++++++++++++++++++++++++++++++++++++++++++++++++ requirements.txt | 4 ---- 2 files changed, 51 insertions(+), 4 deletions(-) create mode 100644 pyproject.toml delete mode 100644 requirements.txt diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..aafe508 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,51 @@ +[build-system] +requires = ["setuptools>=61.0"] +build-backend = "setuptools.build_meta" + +[project] +name = "machine-design" +version = "0.0.1" +authors = [ + { name="Lukáš Adam", email="lukas.adam.cr@gmail.com" }, +] +maintainers = [ + { name="Lukáš Adam", email="lukas.adam.cr@gmail.com" }, +] +description = "Library for designing barriers in electrical machines" +readme = "README.md" +requires-python = ">=3.10" +license = "MIT" +keywords = [ + "design", + "optimisation", +] +classifiers = [ + "Programming Language :: Python :: 3", + "Operating System :: OS Independent", +] +dependencies=[ + 'pyaedt==0.11.2', + 'pyvista[jupyter]', + 'trame', + 'matplotlib', + 'botorch', +] + +[tool.setuptools] +packages = [ + "machine_design", +] + +[tool.ruff] +line-length = 120 + +[tool.ruff.lint] +select = ["E", "F", "I", "UP"] +ignore = ["E501", "F403", "F405"] + +[tool.ruff.format] +quote-style = "double" +indent-style = "space" + +[tool.ruff.lint.per-file-ignores] +"__init__.py" = ["F401"] \ No newline at end of file diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 7c245b5..0000000 --- a/requirements.txt +++ /dev/null @@ -1,4 +0,0 @@ -pyaedt==0.11.2 -pyvista[jupyter] -trame -matplotlib \ No newline at end of file From 2e9cd30ea8db968cb96e8fd30204a36dd0ee63fe Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:54:33 +0000 Subject: [PATCH 52/55] Refactor --- run.py => notebooks/run.py | 0 run_optimization.py => notebooks/run_optimization.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename run.py => notebooks/run.py (100%) rename run_optimization.py => notebooks/run_optimization.py (100%) diff --git a/run.py b/notebooks/run.py similarity index 100% rename from run.py rename to notebooks/run.py diff --git a/run_optimization.py b/notebooks/run_optimization.py similarity index 100% rename from run_optimization.py rename to notebooks/run_optimization.py From 177d4b5719f8dd6fc3b50f1d44159b1c41a3f562 Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:57:06 +0000 Subject: [PATCH 53/55] Ruff format --- machine_design/__init__.py | 10 ++- machine_design/design.py | 139 ++++++++++++++++------------------ machine_design/generators.py | 127 ++++++++++++++++--------------- machine_design/geometry.py | 17 +++-- notebooks/run.py | 46 ++++++----- notebooks/run_optimization.py | 70 +++++++++-------- 6 files changed, 218 insertions(+), 191 deletions(-) diff --git a/machine_design/__init__.py b/machine_design/__init__.py index 3839c2c..abd45fb 100644 --- a/machine_design/__init__.py +++ b/machine_design/__init__.py @@ -1,3 +1,9 @@ from .design import Design -from .generators import HacklGenerator_OneLambda, HacklGenerator_OneLambdaTheta, HacklGenerator_TwoLambdas, FourStupid, save_params -from .geometry import analyze_results, plot_barriers, rotate \ No newline at end of file +from .generators import ( + FourStupid, + HacklGenerator_OneLambda, + HacklGenerator_OneLambdaTheta, + HacklGenerator_TwoLambdas, + save_params, +) +from .geometry import analyze_results, plot_barriers, rotate diff --git a/machine_design/design.py b/machine_design/design.py index 819e974..8234366 100644 --- a/machine_design/design.py +++ b/machine_design/design.py @@ -1,8 +1,8 @@ -import os import numpy as np from ansys.aedt.core import Desktop, Maxwell2d from ansys.aedt.core.modeler.modeler_2d import Modeler2D + class Design: def __init__(self, m2d: Maxwell2d) -> None: self.set_parameters() @@ -10,12 +10,7 @@ def __init__(self, m2d: Maxwell2d) -> None: @classmethod def create(cls, project_name: str, design_name: str, file_name: str, **kwargs) -> "Design": - m2d = Maxwell2d( - project=project_name, - design=design_name, - solution_type="TransientXY", - **kwargs - ) + m2d = Maxwell2d(project=project_name, design=design_name, solution_type="TransientXY", **kwargs) obj = cls(m2d) obj.create_stator() obj.save_project(file_name) @@ -30,9 +25,9 @@ def load(cls, file_name: str, **kwargs) -> "Design": m2d.set_active_design("Design01") return cls(m2d) - + def set_parameters(self) -> None: - #materials + # materials self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" # main definitions @@ -43,7 +38,7 @@ def set_parameters(self) -> None: "SlotNumber": "36", "SlotType": "3", "DiaShaft": "25mm", - "StackLength": "85mm" + "StackLength": "85mm", } # stator slot self.slot_params = { @@ -63,40 +58,47 @@ def set_parameters(self) -> None: "CoilPitch": "9", # coil pitch in slots "SlotLiner": "0.3mm", "SpaceLayers": "0.2mm", - "Nc": "68" # turns per coil + "Nc": "68", # turns per coil } # model parameters PolePairs = 2 - f = 50 #[Hz] - RotSpeed = 60*f/PolePairs #[rpm] + f = 50 # [Hz] + RotSpeed = 60 * f / PolePairs # [rpm] self.mod_params = { "Poles": f"2*{PolePairs}", "ModelLength": "85mm", "SymmetryFactor": "Poles", "StatorSkewAngle": "0deg", } - #operation parameters + # operation parameters self.oper_params = { "Im": "1.5*sqrt(2)A", - "epsI": "pi/4", #current angle + "epsI": "pi/4", # current angle "InitPos": "-30deg", "f": f"{f}Hz", "RotSpeed": f"{RotSpeed}rpm", - "Nper": "1/6" - "", #number of included periods - "PointPer": "101" #number of time points per period + "Nper": "1/6", # number of included periods + "PointPer": "101", # number of time points per period } self.rot_points = [ ["DiaShaft/2*cos(360deg/SymmetryFactor)", "DiaShaft/2*sin(360deg/SymmetryFactor)", "0mm"], ["DiaShaft/2*cos(360deg/(2*SymmetryFactor))", "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", "0mm"], ["DiaShaft/2", "0mm", "0mm"], ["DiaStatorGap/2-Airgap", "0mm", "0mm"], - ["(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", "0mm"], - ["(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", "0mm"], + [ + "(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", + "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", + "0mm", + ], + [ + "(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", + "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", + "0mm", + ], ] self.setup_name = "Setup1" - self.rotor_r_min = self.mm_to_str('geom_params', 'DiaShaft')/2 - self.rotor_r_max = self.mm_to_str('geom_params', 'DiaStatorGap')/2 - self.mm_to_str('geom_params', 'Airgap') + self.rotor_r_min = self.mm_to_str("geom_params", "DiaShaft") / 2 + self.rotor_r_max = self.mm_to_str("geom_params", "DiaStatorGap") / 2 - self.mm_to_str("geom_params", "Airgap") def create_stator(self) -> None: m2d = self.m2d @@ -141,7 +143,7 @@ def create_stator(self) -> None: name="Shaft", ) - #motion setup + # motion setup m2d.assign_rotate_motion( assignment="Band", coordinate_system="Global", @@ -199,7 +201,7 @@ def create_stator(self) -> None: stator_id.color = (192, 192, 192) # rgb stator_id.transparency = 0.0 - #Winding + # Winding coil_id = modeler.create_rectangle( origin=["DiaStatorGap/2+Hs0+Hs1+SlotLiner", "-(Bs1/2-SlotLiner)", 0], sizes=["(Hs2+Rs/2-2*SlotLiner)/Layers-(Layers-1)*SpaceLayers/2", "Bs1-2*SlotLiner", 0], @@ -209,9 +211,7 @@ def create_stator(self) -> None: coil_id.color = (255, 128, 0) coil_id.transparency = 0.0 modeler.rotate(assignment=coil_id, axis="Z", angle="360deg/SlotNumber/2") - coil_id.duplicate_around_axis( - axis="Z", angle="360deg/SlotNumber", clones="CoilPitch", create_new_objects=True - ) + coil_id.duplicate_around_axis(axis="Z", angle="360deg/SlotNumber", clones="CoilPitch", create_new_objects=True) id_coils = modeler.get_objects_w_string(string_name="Coil", case_sensitive=True) # Create section of machine @@ -231,9 +231,7 @@ def create_stator(self) -> None: # Symmetrical boundary conditions pos_1 = "DiaStatorGap/4" - id_bc_1 = modeler.get_edgeid_from_position( - position=[pos_1, 0, 0], assignment="Region" - ) + id_bc_1 = modeler.get_edgeid_from_position(position=[pos_1, 0, 0], assignment="Region") id_bc_2 = modeler.get_edgeid_from_position( position=[ pos_1 + "*cos((360deg/SymmetryFactor))", @@ -260,16 +258,14 @@ def create_stator(self) -> None: ], assignment="Region", ) - m2d.assign_vector_potential( - assignment=id_bc_az, vector_value=0, boundary="A0" - ) + m2d.assign_vector_potential(assignment=id_bc_az, vector_value=0, boundary="A0") # Excitations I_A = "Im * cos(2*pi*f*time+epsI)" I_B = "Im * cos(2*pi*f*time-120deg+epsI)" I_C = "Im * cos(2*pi*f*time-240deg+epsI)" - #Define phase windings + # Define phase windings m2d.assign_coil( assignment=["Coil"], conductors_number="Nc", @@ -296,9 +292,7 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseA", ) - m2d.add_winding_coils( - assignment="PhaseA", coils=["CS1", "CS2", "CS3"] - ) + m2d.add_winding_coils(assignment="PhaseA", coils=["CS1", "CS2", "CS3"]) m2d.assign_coil( assignment=["Coil_6"], conductors_number="Nc", @@ -325,9 +319,7 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseB", ) - m2d.add_winding_coils( - assignment="PhaseB", coils=["CS7", "CS8", "CS9"] - ) + m2d.add_winding_coils(assignment="PhaseB", coils=["CS7", "CS8", "CS9"]) m2d.assign_coil( assignment=["Coil_3"], conductors_number="Nc", @@ -354,9 +346,7 @@ def create_stator(self) -> None: parallel_branches="ParallelPaths", name="PhaseC", ) - m2d.add_winding_coils( - assignment="PhaseC", coils=["CS4", "CS5", "CS6"] - ) + m2d.add_winding_coils(assignment="PhaseC", coils=["CS4", "CS5", "CS6"]) # Mesh operation m2d.mesh.assign_length_mesh( @@ -374,18 +364,16 @@ def create_stator(self) -> None: name="stator", ) - #core loss + # core loss m2d.set_core_losses("Stator", core_loss_on_field=False) - #inductance calculation - m2d.change_inductance_computation( - compute_transient_inductance=True, incremental_matrix=False - ) - #model depth + # inductance calculation + m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=False) + # model depth m2d.model_depth = "StackLength" - #symmetry + # symmetry m2d.change_symmetry_multiplier("SymmetryFactor") - # Calculation setup + # Calculation setup setup = m2d.create_setup(name=self.setup_name) setup.props["StopTime"] = "Nper/f" setup.props["TimeStep"] = "1/(f*(PointPer-1))" @@ -396,7 +384,7 @@ def create_stator(self) -> None: setup.update() m2d.validate_simple() - #ooutput variables + # ooutput variables output_vars = { "pos": "(Moving1.Position -InitPos) * Poles/2", "cos0": "cos(pos)", @@ -424,22 +412,25 @@ def create_stator(self) -> None: for k, v in output_vars.items(): m2d.create_output_variable(k, v) - #Definitions for plots + # Definitions for plots post_params = { # reports - ("InducedVoltage(PhaseA)","InducedVoltage(PhaseB)","InducedVoltage(PhaseC)"): "InducedVoltage", + ("InducedVoltage(PhaseA)", "InducedVoltage(PhaseB)", "InducedVoltage(PhaseC)"): "InducedVoltage", ("Moving1.Torque"): "Torque", - ("InputCurrent(PhaseA)","InputCurrent(PhaseB)","InputCurrent(PhaseC)"): "Current", - ("FluxLinkage(PhaseA)","FluxLinkage(PhaseB)","FluxLinkage(PhaseC)",): "FluxLinkage", + ("InputCurrent(PhaseA)", "InputCurrent(PhaseB)", "InputCurrent(PhaseC)"): "Current", + ( + "FluxLinkage(PhaseA)", + "FluxLinkage(PhaseB)", + "FluxLinkage(PhaseC)", + ): "FluxLinkage", ("I_d", "I_q"): "Current_dq", ("Flux_d", "Flux_q"): "FluxLinkage_dq", ("Ui_d", "Ui_q"): "InducedVoltage_dq", ("L_d", "L_q"): "Inductance_dq", - } - #Create Report + # Create Report for k, v in post_params.items(): - expressions = list(k) if isinstance(k, tuple) else [k] #if multiple report, use list(k). Else, use k - report = m2d.post.create_report( + expressions = list(k) if isinstance(k, tuple) else [k] # if multiple report, use list(k). Else, use k + m2d.post.create_report( expressions=expressions, setup_sweep_name="", domain="Sweep", @@ -459,7 +450,7 @@ def add_rotor(self) -> None: assert isinstance(modeler, Modeler2D) rotor_id = modeler.create_polyline( - points=self.rot_points, segment_type=[ "Arc","Line", "Arc"], cover_surface=True, name="Rotor" + points=self.rot_points, segment_type=["Arc", "Line", "Arc"], cover_surface=True, name="Rotor" ) self.rotor_id = rotor_id rotor_id.material_name = self.Fe @@ -473,10 +464,10 @@ def add_rotor_barrier(self, barrier_points, segment_type=None) -> None: # Round it for Ansys barrier_points = np.round(barrier_points, 6) - # Potentially add the z axis + # Potentially add the z axis if barrier_points.shape[1] == 2: - barrier_points = np.hstack((barrier_points, np.zeros((len(barrier_points),1)))) - + barrier_points = np.hstack((barrier_points, np.zeros((len(barrier_points), 1)))) + # Convert them into a string format and interpolate points_str = [[str(y) for y in x] for x in barrier_points] barrier_id = modeler.create_polyline( @@ -499,27 +490,25 @@ def compute(self, NUM_CORES: int = 1): maximum_elements=None, name="rotor", ) - #core loss rotor + # core loss rotor m2d.set_core_losses("Rotor", core_loss_on_field=False) # Analyze m2d.analyze_setup(self.setup_name, use_auto_settings=False, cores=NUM_CORES) - solutions = m2d.post.get_solution_data( - expressions="Moving1.Torque", primary_sweep_variable="Time" - ) + solutions = m2d.post.get_solution_data(expressions="Moving1.Torque", primary_sweep_variable="Time") try: return solutions.data_magnitude() except AttributeError: return None - + def delete_rotor(self) -> None: assert isinstance(self.m2d.modeler, Modeler2D) self.m2d.modeler.delete(self.rotor_id) def save_design(self, file_name: str, **kwargs) -> None: - show = kwargs.pop('show', False) - view = kwargs.pop('view', 'xy') + show = kwargs.pop("show", False) + view = kwargs.pop("view", "xy") self.m2d.plot(show=show, output_file=file_name, view=view) def save_project(self, file_name: str | None = None) -> None: @@ -527,12 +516,12 @@ def save_project(self, file_name: str | None = None) -> None: self.m2d.save_project() else: self.m2d.save_project(file_name) - + def close_project(self) -> None: self.m2d.close_desktop() - + def mm_to_str(self, var, field) -> float: val = getattr(self, var)[field] - if not val.endswith('mm'): - raise Exception('val must end with mm') + if not val.endswith("mm"): + raise Exception("val must end with mm") return float(val[:-2]) diff --git a/machine_design/generators.py b/machine_design/generators.py index f24cb32..fd271bd 100644 --- a/machine_design/generators.py +++ b/machine_design/generators.py @@ -1,11 +1,12 @@ +import pickle from abc import ABC, abstractmethod +from typing import Any + import numpy as np -from scipy.interpolate import CubicSpline -import pickle -from shapely.geometry import LineString from design import Design from geometry import rotate -from typing import Any +from scipy.interpolate import CubicSpline +from shapely.geometry import LineString def get_arc(R: float, start_deg: float, end_deg: float, n_points: int) -> np.ndarray: @@ -16,28 +17,33 @@ def get_arc(R: float, start_deg: float, end_deg: float, n_points: int) -> np.nda return np.column_stack((x, y)) + def signed_distance(P: np.ndarray, a: float, b: float, c: float) -> np.ndarray: - return a * P[:,0] + b * P[:,1] + c + return a * P[:, 0] + b * P[:, 1] + c + def check_barrier(barrier: np.ndarray) -> None: if not np.array_equal(barrier[0], barrier[-1]): raise ValueError("Barrier must start and end at the same point") + def save_params(params, file_name: str) -> None: - with open(file_name, 'wb') as handle: + with open(file_name, "wb") as handle: pickle.dump(params, handle, protocol=pickle.HIGHEST_PROTOCOL) class BarrierGenerator(ABC): # TODO: rename offset - def __init__(self, design: Design, r_stator_end: float, offset: None | float = None, n_curve: int = 500, n_flat: int = 20) -> None: + def __init__( + self, design: Design, r_stator_end: float, offset: None | float = None, n_curve: int = 500, n_flat: int = 20 + ) -> None: self.r_max = design.rotor_r_max self.r_min = design.rotor_r_min self.R = self.r_max - r_stator_end self.offset = offset self.n_curve = n_curve self.n_flat = n_flat - + @abstractmethod def generate_barriers(self) -> list[np.ndarray]: pass @@ -72,10 +78,10 @@ def feasible_barriers(self, barriers: list[np.ndarray], tol: float = 1e-8) -> bo for curve in curves: if not curve.is_simple: return False - + # Check if two curves intersect themselves for i in range(n): - for j in range(i+1,n): + for j in range(i + 1, n): if curves[i].intersects(curves[j]): return False return True @@ -86,7 +92,7 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: lam = np.linspace(0, 1, self.n_flat)[:, None] barriers = [] - for a, b, c in zip([1., -1.], [-1., 1.], [-self.offset, -self.offset]): + for a, b, c in zip([1.0, -1.0], [-1.0, 1.0], [-self.offset, -self.offset]): f = signed_distance(barrier, a, b, c) f_is_positive = f > 0 mask_positive = np.where(f_is_positive)[0] @@ -94,7 +100,7 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: barrier_new = [barrier[mask_positive[0]]] for i in range(len(mask_positive) - 1): i1 = mask_positive[i] - i2 = mask_positive[i+1] + i2 = mask_positive[i + 1] if i2 - i1 == 1: barrier_new.append(barrier[i2]) @@ -114,8 +120,9 @@ def split_barrier(self, barrier: np.ndarray) -> list[np.ndarray]: def split_barriers(self, barriers: list[np.ndarray]) -> list[np.ndarray]: return [x for barrier in barriers for x in self.split_barrier(barrier)] + class FourStupid(BarrierGenerator): - def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwargs) -> None: + def __init__(self, design, r_stator_end, der1=1.0, der2=1.0, symmetric=True, **kwargs) -> None: self.der1 = der1 self.der2 = der2 self.w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0 @@ -123,31 +130,31 @@ def __init__(self, design, r_stator_end, der1=1., der2=1., symmetric=True, **kwa super().__init__(design, r_stator_end, **kwargs) def _create_barrier( - self, - y_min, - w_min, - y_mid, - w_mid, - theta, - w_max, - ): - - theta1 = (theta+45) / 180*np.pi - x_max1 = self.R*np.cos(theta1) - y_max1 = self.R*np.sin(theta1) - - x1 = [0, x_max1/2, x_max1] + self, + y_min, + w_min, + y_mid, + w_mid, + theta, + w_max, + ): + + theta1 = (theta + 45) / 180 * np.pi + x_max1 = self.R * np.cos(theta1) + y_max1 = self.R * np.sin(theta1) + + x1 = [0, x_max1 / 2, x_max1] y1 = [y_min, y_mid, y_max1] - f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1,self.der1))) + f1 = CubicSpline(x1, y1, bc_type=((1, 0), (1, self.der1))) theta2 = theta1 + w_max / self.r_max - x_max2 = self.R*np.cos(theta2) - y_max2 = self.R*np.sin(theta2) + x_max2 = self.R * np.cos(theta2) + y_max2 = self.R * np.sin(theta2) - s = w_mid / np.sqrt(1 + self.der1**2/4) - x2 = [0, x1[1]-s*self.der1/2, x_max2] - y2 = [y_min+w_min, y1[1]+s, y_max2] - f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1,self.der2))) + s = w_mid / np.sqrt(1 + self.der1**2 / 4) + x2 = [0, x1[1] - s * self.der1 / 2, x_max2] + y2 = [y_min + w_min, y1[1] + s, y_max2] + f2 = CubicSpline(x2, y2, bc_type=((1, 0), (1, self.der2))) x_interp1 = np.linspace(x1[0], x1[-1], self.n_curve) x_interp2 = np.linspace(x2[0], x2[-1], self.n_curve) @@ -167,20 +174,20 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: return lb, ub def random_parameters(self): - rand1 = 0.5*np.random.random(4) - rand2 = 13+4*np.random.random() - rand3 = 2+3*np.random.random() - rand4 = 2+3*np.random.random() - rand5 = 2+2*np.random.random() + rand1 = 0.5 * np.random.random(4) + rand2 = 13 + 4 * np.random.random() + rand3 = 2 + 3 * np.random.random() + rand4 = 2 + 3 * np.random.random() + rand5 = 2 + 2 * np.random.random() return rand1, rand2, rand3, rand4, rand5 def set_parameters(self, params) -> None: rand1, rand2, rand3, rand4, rand5 = params self.w_mins = self.w_mins_base + rand1 y_min0 = rand2 - y_min1 = y_min0+self.w_mins[0] + rand3 - y_min2 = y_min1+self.w_mins[1] + rand4 - y_min3 = y_min2+self.w_mins[2] + rand5 + y_min1 = y_min0 + self.w_mins[0] + rand3 + y_min2 = y_min1 + self.w_mins[1] + rand4 + y_min3 = y_min2 + self.w_mins[2] + rand5 self.w_mids = self.w_mins self.y_mins = np.array([y_min0, y_min1, y_min2, y_min3]) self.y_mids = self.y_mins + np.array([2, 1.5, 1, 0.5]) @@ -197,14 +204,14 @@ def generate_barriers(self) -> list[np.ndarray]: xy_all = np.vstack((x_all, y_all)).T barriers.append(xy_all) return barriers - + class AbstractHacklGenerator(BarrierGenerator): def __init__(self, design, r_stator_end, **kwargs): - self.phis_inner_min = np.asarray([4., 17.3, 30.6]) + self.phis_inner_min = np.asarray([4.0, 17.3, 30.6]) self.phis_inner_max = np.asarray([10.2, 23.5, 36.8]) - self.phis_outer_min = np.asarray([10.6, 24., 37.3]) - self.phis_outer_max = np.asarray([16.8, 30.2, 44.]) + self.phis_outer_min = np.asarray([10.6, 24.0, 37.3]) + self.phis_outer_max = np.asarray([16.8, 30.2, 44.0]) self.n_barriers = len(self.phis_inner_min) super().__init__(design, r_stator_end, **kwargs) @@ -213,8 +220,8 @@ def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple pass def random_parameters(self): - phis_inner = self.phis_inner_min + (self.phis_inner_max-self.phis_inner_min)*np.random.rand(self.n_barriers) - phis_outer = self.phis_outer_min + (self.phis_outer_max-self.phis_outer_min)*np.random.rand(self.n_barriers) + phis_inner = self.phis_inner_min + (self.phis_inner_max - self.phis_inner_min) * np.random.rand(self.n_barriers) + phis_outer = self.phis_outer_min + (self.phis_outer_max - self.phis_outer_min) * np.random.rand(self.n_barriers) return phis_inner, phis_outer def set_parameters(self, params) -> None: @@ -228,11 +235,11 @@ def generate_barriers(self) -> list[np.ndarray]: # Generate the longer part of the barrier pts_outer = self._get_bezier_curve(phi_outer, False, i) pts_inner = self._get_bezier_curve(phi_inner, True, i) - + # Generate the arcs at the end arc_top = get_arc(self.R, 90 - phi_outer, 90 - phi_inner, self.n_flat) arc_bottom = get_arc(self.R, phi_inner, phi_outer, self.n_flat) - + # Merge them together barrier = np.concatenate((pts_outer, arc_top[1:-1], pts_inner[::-1], arc_bottom[1:])) barriers.append(barrier) @@ -252,8 +259,8 @@ def _get_bezier_curve(self, phi_deg, is_inner, i): # Formula 33: Polynomial expansion to generate continuous points z_vals = np.linspace(0, 1, self.n_curve)[:, None] return ( - (1 - z_vals)**3 * r0 - + 3 * (1 - z_vals)**2 * z_vals * r1 + (1 - z_vals) ** 3 * r0 + + 3 * (1 - z_vals) ** 2 * z_vals * r1 + 3 * (1 - z_vals) * z_vals**2 * r2 + z_vals**3 * r3 ) @@ -270,12 +277,12 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: lb = np.concatenate((self.phis_inner_min, self.phis_outer_min, [self.lam_min])) ub = np.concatenate((self.phis_inner_max, self.phis_outer_max, [self.lam_max])) return lb, ub - + def random_parameters(self): pars = super().random_parameters() - lam = self.lam_min + (self.lam_max-self.lam_min)*np.random.rand(1)[0] + lam = self.lam_min + (self.lam_max - self.lam_min) * np.random.rand(1)[0] return *pars, lam - + def set_parameters(self, params) -> None: phis_inner, phis_outer, lam = params super().set_parameters((phis_inner, phis_outer)) @@ -283,7 +290,7 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6] - + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: return self.lam * x + (1 - self.lam) * y, y @@ -302,7 +309,7 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: def random_parameters(self): pars = super().random_parameters() - theta = self.theta_min + (self.theta_max-self.theta_min)*np.random.rand(1)[0] + theta = self.theta_min + (self.theta_max - self.theta_min) * np.random.rand(1)[0] return *pars, theta def set_parameters(self, params) -> None: @@ -312,7 +319,7 @@ def set_parameters(self, params) -> None: def X_to_params(self, X: np.ndarray): return X[:3], X[3:6], X[6], X[7] - + def _bezier_point(self, x: float, y: float, is_inner: bool, order: int) -> tuple[float, float]: x, y = super()._bezier_point(x, y, is_inner, order) if order < 2: @@ -337,8 +344,8 @@ def bounds(self) -> tuple[np.ndarray, np.ndarray]: def random_parameters(self): pars = super().random_parameters() - lam_inner = self.lam_inner_min + (self.lam_inner_max-self.lam_inner_min)*np.random.rand(1)[0] - lam_outer = self.lam_outer_min + (self.lam_outer_max-self.lam_outer_min)*np.random.rand(1)[0] + lam_inner = self.lam_inner_min + (self.lam_inner_max - self.lam_inner_min) * np.random.rand(1)[0] + lam_outer = self.lam_outer_min + (self.lam_outer_max - self.lam_outer_min) * np.random.rand(1)[0] return *pars, lam_inner, lam_outer def set_parameters(self, params) -> None: diff --git a/machine_design/geometry.py b/machine_design/geometry.py index b67ec07..ecac422 100644 --- a/machine_design/geometry.py +++ b/machine_design/geometry.py @@ -1,24 +1,26 @@ -import numpy as np import matplotlib.pyplot as plt +import numpy as np def rotation_matrix(alpha, rad=False): if not rad: - alpha /= 360/(2*np.pi) + alpha /= 360 / (2 * np.pi) return np.array([[np.cos(alpha), -np.sin(alpha)], [np.sin(alpha), np.cos(alpha)]]) + def rotate(x, y, alpha, **kwargs): T = rotation_matrix(alpha, **kwargs) xy_rot = T @ np.vstack((x, y)) return xy_rot[0], xy_rot[1] + def plot_barriers(barriers, design, title=None, file_name=None): plt.figure() - alphas = np.linspace(0*np.pi,2*np.pi/4,100) + alphas = np.linspace(0 * np.pi, 2 * np.pi / 4, 100) for r in [design.rotor_r_min, design.rotor_r_max]: - plt.plot(r*np.cos(alphas), r*np.sin(alphas)) + plt.plot(r * np.cos(alphas), r * np.sin(alphas)) for barrier in barriers: - plt.plot(barrier[:,0], barrier[:,1]) + plt.plot(barrier[:, 0], barrier[:, 1]) plt.axis("equal") if title is not None: plt.title(title) @@ -26,8 +28,9 @@ def plot_barriers(barriers, design, title=None, file_name=None): plt.savefig(file_name) plt.close() + def analyze_results(Tor: np.ndarray) -> tuple[float, float, float]: TorAvg = np.mean(Tor[:-1]) - TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1]-TorAvg))) - TorRippleRms = TorRmsAC/TorAvg*100 + TorRmsAC = np.sqrt(np.mean(np.square(Tor[:-1] - TorAvg))) + TorRippleRms = TorRmsAC / TorAvg * 100 return TorAvg, TorRmsAC, TorRippleRms diff --git a/notebooks/run.py b/notebooks/run.py index 1ca3622..b19e2b1 100644 --- a/notebooks/run.py +++ b/notebooks/run.py @@ -1,30 +1,42 @@ # export ANSYSEM_ROOT241=/data/AnsysEM/v241/Linux64 import os + import numpy as np import pandas as pd -from machine_design import Design, analyze_results, plot_barriers -from machine_design import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta, save_params + +from machine_design import ( + Design, + FourStupid, + HacklGenerator_OneLambda, + HacklGenerator_OneLambdaTheta, + HacklGenerator_TwoLambdas, + analyze_results, + plot_barriers, + save_params, +) project_name = "SynRM_test" design_name = "Design01" -path_data = os.path.join(os.getcwd(), 'data') -path_results = 'results' +path_data = os.path.join(os.getcwd(), "data") +path_results = "results" for path in [path_data, path_results]: os.makedirs(path, exist_ok=True) -file_name_aedt = f'{path_data}/{project_name}.aedt' +file_name_aedt = f"{path_data}/{project_name}.aedt" plot_design = True n_designs = 20 # Define constants AEDT_VERSION = "2024.1" NUM_CORES = 4 -NG_MODE = True #non-graphical mode -CLS_EXIT = True #close on exit +NG_MODE = True # non-graphical mode +CLS_EXIT = True # close on exit if not os.path.exists(file_name_aedt): design = Design.create( - project_name, design_name, file_name_aedt, + project_name, + design_name, + file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, new_desktop=False, @@ -64,7 +76,7 @@ Tor = design.compute(NUM_CORES) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan - else: + else: TorAvg, _, TorRippleRms = analyze_results(Tor) # Delete the rotor @@ -72,18 +84,18 @@ # Potentially save the design if plot_design: - title = f'Torque mean value: {np.round(TorAvg,2)} Nm, ripple relative value: {np.round(TorRippleRms,2)} %' - file_name = f'{path_results}/design_{generator.name}_{i}' - plot_barriers(barriers, design, title=title, file_name = f"{file_name}.png") + title = f"Torque mean value: {np.round(TorAvg, 2)} Nm, ripple relative value: {np.round(TorRippleRms, 2)} %" + file_name = f"{path_results}/design_{generator.name}_{i}" + plot_barriers(barriers, design, title=title, file_name=f"{file_name}.png") save_params(params, f"{file_name}.pkl") metadata_new = { - 'method': generator.name, - 'design': i, - 'T': TorAvg, - 'ripple': TorRippleRms, + "method": generator.name, + "design": i, + "T": TorAvg, + "ripple": TorRippleRms, } metadata = pd.concat((metadata, pd.DataFrame([metadata_new])), ignore_index=True) metadata.to_csv(f"{path_results}/metadata.csv", index=False) -design.close_project() \ No newline at end of file +design.close_project() diff --git a/notebooks/run_optimization.py b/notebooks/run_optimization.py index 21b22b3..d7e28b2 100644 --- a/notebooks/run_optimization.py +++ b/notebooks/run_optimization.py @@ -1,37 +1,43 @@ # export ANSYSEM_ROOT241=/data/AnsysEM/v241/Linux64 -from collections.abc import Iterable import os +import pickle +from collections.abc import Iterable + import numpy as np import pandas as pd -import pickle -from machine_design import Design, analyze_results, plot_barriers -from machine_design import FourStupid, HacklGenerator_OneLambda, HacklGenerator_TwoLambdas, HacklGenerator_OneLambdaTheta -import torch -torch.set_default_dtype(torch.float64) -from torch import Tensor -from botorch.utils.transforms import normalize, unnormalize import torch - -from botorch.models import SingleTaskGP from botorch import fit_gpytorch_mll -from gpytorch.mlls import ExactMarginalLogLikelihood - from botorch.acquisition.multi_objective.logei import qLogExpectedHypervolumeImprovement +from botorch.models import SingleTaskGP +from botorch.optim import optimize_acqf from botorch.utils.multi_objective.box_decompositions import NondominatedPartitioning from botorch.utils.multi_objective.pareto import is_non_dominated -from botorch.optim import optimize_acqf +from botorch.utils.transforms import normalize, unnormalize +from gpytorch.mlls import ExactMarginalLogLikelihood +from torch import Tensor + +from machine_design import ( + Design, + HacklGenerator_OneLambda, + analyze_results, +) + +torch.set_default_dtype(torch.float64) + def objective_transform(TorAvg, TorRippleRms): if pd.isnull(TorAvg): return -99999, -99999 else: - return TorAvg, -TorRippleRms/100 + return TorAvg, -TorRippleRms / 100 + def objective(Xs: Tensor, *args) -> Tensor: vals = [objective_single(X, *args) for X in Xs] return torch.stack(vals, dim=0) + def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: X = unnormalize(X, bounds) params = generator.X_to_params(X.numpy()) @@ -49,7 +55,7 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: Tor = design.compute(NUM_CORES) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan - else: + else: TorAvg, _, TorRippleRms = analyze_results(Tor) # Delete the rotor @@ -58,6 +64,7 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES) -> Tensor: f1, f2 = objective_transform(TorAvg, TorRippleRms) return torch.tensor([f1, f2]) + def init_points(root, method): results = pd.read_csv(f"{root}/metadata.csv") results = results[~results["T"].isnull()] @@ -66,9 +73,9 @@ def init_points(root, method): Xs, Ys = [], [] for _, r in results.iterrows(): - with open(r["path"], 'rb') as f: + with open(r["path"], "rb") as f: params = pickle.load(f) - + X = [] for x in params: if isinstance(x, Iterable): @@ -83,21 +90,24 @@ def init_points(root, method): return torch.stack(Xs), torch.stack(Ys) + project_name = "SynRM_test" design_name = "Design01" -path_data = os.path.join(os.getcwd(), 'data') +path_data = os.path.join(os.getcwd(), "data") os.makedirs(path_data, exist_ok=True) -file_name_aedt = f'{path_data}/{project_name}.aedt' +file_name_aedt = f"{path_data}/{project_name}.aedt" # Define constants AEDT_VERSION = "2024.1" NUM_CORES = 4 -NG_MODE = True #non-graphical mode -CLS_EXIT = True #close on exit +NG_MODE = True # non-graphical mode +CLS_EXIT = True # close on exit if not os.path.exists(file_name_aedt): design = Design.create( - project_name, design_name, file_name_aedt, + project_name, + design_name, + file_name_aedt, version=AEDT_VERSION, non_graphical=NG_MODE, new_desktop=False, @@ -127,15 +137,18 @@ def init_points(root, method): train_X, train_Y = init_points(root_init, method) train_X = normalize(train_X, bounds) bounds_normalized = normalize(bounds, bounds) -ref_point = torch.tensor([3.8,-max_ripple]) +ref_point = torch.tensor([3.8, -max_ripple]) + def objective_lambda(Xs): return objective(Xs, design, generator, bounds, NUM_CORES) + def ripple_constraint(Y): - ripple = -Y[...,1] + ripple = -Y[..., 1] return ripple - max_ripple + for _ in range(n_iters): # Fit surrogate model = SingleTaskGP(train_X, train_Y) @@ -144,17 +157,14 @@ def ripple_constraint(Y): # Compute Pareto front pareto_Y = train_Y[is_non_dominated(train_Y)] - partitioning = NondominatedPartitioning( - ref_point=ref_point, - Y=pareto_Y - ) + partitioning = NondominatedPartitioning(ref_point=ref_point, Y=pareto_Y) # Define acquisition function acq = qLogExpectedHypervolumeImprovement( model=model, ref_point=ref_point.tolist(), partitioning=partitioning, - constraints=[ripple_constraint], + constraints=[ripple_constraint], ) # Optimize acquisition function to select candidate points. Reject unfeasible points @@ -189,4 +199,4 @@ def ripple_constraint(Y): train_Y = torch.cat([train_Y, new_Y]) np.savez(f"results_{method}.npz", train_X=unnormalize(train_X, bounds), train_Y=train_Y) -design.close_project() \ No newline at end of file +design.close_project() From 34f8ffd61f510d3ad62b33d4c037d5b85bc3ce4b Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:57:37 +0000 Subject: [PATCH 54/55] Notebook added --- notebooks/analyze.ipynb | 177 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 notebooks/analyze.ipynb diff --git a/notebooks/analyze.ipynb b/notebooks/analyze.ipynb new file mode 100644 index 0000000..9e18ba7 --- /dev/null +++ b/notebooks/analyze.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "652452a4", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import torch\n", + "from botorch.utils.multi_objective import Hypervolume\n", + "from botorch.utils.multi_objective.pareto import is_non_dominated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb2df116", + "metadata": {}, + "outputs": [], + "source": [ + "def get_pareto(Y):\n", + " idx_pareto = is_non_dominated(torch.from_numpy(Y))\n", + " idx_sort = np.argsort(Y[idx_pareto, 0])\n", + " idx_combined = np.where(idx_pareto)[0][idx_sort]\n", + "\n", + " return Y[idx_combined], idx_combined\n", + "\n", + "\n", + "def load_data(file_name):\n", + " data = np.load(file_name)\n", + " X = data[\"train_X\"]\n", + " Y = data[\"train_Y\"]\n", + " pareto_Y, idx_pareto = get_pareto(Y)\n", + " Y[:, 1] = -100 * Y[:, 1]\n", + " pareto_Y[:, 1] = -100 * pareto_Y[:, 1]\n", + "\n", + " return X, Y, pareto_Y, idx_pareto" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4bf82b18", + "metadata": {}, + "outputs": [], + "source": [ + "n_init = 20\n", + "colors = [\"black\", \"blue\", \"orange\", \"green\"]\n", + "methods_rename = {\n", + " \"HacklGenerator_OneLambda\": (\"Same lambda\", \"black\"),\n", + " \"HacklGenerator_TwoLambdas\": (\"Different lambda\", \"blue\"),\n", + " \"FourStupid\": (\"Proposed method\", \"orange\"),\n", + "}\n", + "labels = [\"All designs\", \"Initial designs\"]\n", + "methods = list(methods_rename.keys())\n", + "prefixes = [\"\", \"old_\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d6eaacb", + "metadata": {}, + "outputs": [], + "source": [ + "for prefix in prefixes:\n", + " X = {}\n", + " Y = {}\n", + " pareto_X = {}\n", + " pareto_Y = {}\n", + " for prefix in prefixes:\n", + " for method in methods:\n", + " key = f\"{prefix}{method}\"\n", + " file_name = f\"results/{prefix}results_{method}.npz\"\n", + " X[key], Y[key], _, idx = load_data(file_name)\n", + " pareto_X[key] = X[key][idx]\n", + " pareto_Y[key] = Y[key][idx]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2bf03f9", + "metadata": {}, + "outputs": [], + "source": [ + "max_ripple = 0.1\n", + "ref_point = [3.8, -max_ripple]\n", + "hypervolume = Hypervolume(torch.tensor(ref_point))\n", + "\n", + "for i, prefix in enumerate(prefixes):\n", + " Y_all = np.vstack([Y[f\"{prefix}{method}\"] for method in methods])\n", + " Y_init = np.vstack([Y[f\"{prefix}{method}\"][:n_init] for method in methods])\n", + " Y_orig = Y_all.copy()\n", + " Y_orig[:, 1] = -1 / 100 * Y_orig[:, 1]\n", + " _, idx = get_pareto(Y_orig)\n", + " Y_pareto = Y_all[idx]\n", + " hv = hypervolume.compute(torch.from_numpy(Y_orig))\n", + "\n", + " plt.scatter(Y_all[:, 0], Y_all[:, 1], color=\"gray\", alpha=0.2, s=12, zorder=1, label=labels[0] if i == 0 else \"\")\n", + " plt.scatter(Y_init[:, 0], Y_init[:, 1], color=\"blue\", s=14, label=labels[1] if i == 0 else \"\")\n", + " if i == 0:\n", + " color = \"black\"\n", + " label = f\"Pareto front with constraints (HV={100 * hv:.1f})\"\n", + " else:\n", + " color = \"red\"\n", + " label = f\"Pareto front without constraints (HV={100 * hv:.1f})\"\n", + " plt.scatter(Y_pareto[:, 0], Y_pareto[:, 1], color=color, s=14, label=label)\n", + " plt.plot(Y_pareto[:, 0], Y_pareto[:, 1], color=color)\n", + "plt.scatter([ref_point[0]], [-100 * ref_point[1]], color=\"black\", label=\"Reference point\", marker=\"*\")\n", + "\n", + "plt.xlim([3.2, 4.6])\n", + "plt.ylim([0, 20])\n", + "plt.legend()\n", + "plt.grid()\n", + "ax = plt.gca()\n", + "ax.set_axisbelow(True)\n", + "plt.savefig(\"all_designs.png\", dpi=600, bbox_inches=\"tight\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f866acf", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(6, 4))\n", + "for i, prefix in enumerate(prefixes):\n", + " ax = plt.subplot(1, 2, i + 1)\n", + " for method in methods:\n", + " Y_pareto = pareto_Y[f\"{prefix}{method}\"]\n", + " label, color = methods_rename[method]\n", + " plt.plot(Y_pareto[:, 0], Y_pareto[:, 1], label=label, color=color, linewidth=3)\n", + " ax.axhline(10, color=\"black\", linestyle=\":\", linewidth=2, alpha=0.7)\n", + "\n", + " if i == 0:\n", + " plt.legend()\n", + " plt.title(\"With constraints\")\n", + " plt.ylabel(\"RMS ripple [%]\")\n", + " else:\n", + " plt.title(\"Without constraints\")\n", + "\n", + " plt.xlabel(\"Mean torque [Nm]\")\n", + " plt.xlim([4.1, 4.6])\n", + " plt.ylim([2, 20])\n", + " plt.grid(alpha=0.35)\n", + " ax.fill_between([4.1, 4.6], 0, 10, color=\"green\", alpha=0.05)\n", + "plt.savefig(\"result.png\", dpi=600, bbox_inches=\"tight\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ansys", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 3f715b2c4798c6c0d4a054bffd11f534c77af9e3 Mon Sep 17 00:00:00 2001 From: sadda Date: Tue, 17 Mar 2026 09:58:53 +0000 Subject: [PATCH 55/55] Added Ruff workflow --- .github/workflows/rull.yml | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 .github/workflows/rull.yml diff --git a/.github/workflows/rull.yml b/.github/workflows/rull.yml new file mode 100644 index 0000000..3ea32a9 --- /dev/null +++ b/.github/workflows/rull.yml @@ -0,0 +1,24 @@ +name: Ruff + +on: + push: + branches: + - main + pull_request: + +jobs: + ruff: + runs-on: ubuntu-latest + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Install Ruff + run: pip install ruff + + - name: Ruff lint + run: ruff check . + + - name: Ruff format check + run: ruff format . --check \ No newline at end of file