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 diff --git a/.gitignore b/.gitignore index bf911bb..687c547 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,8 @@ data/ -results/ +results*/ +results_*.npz +nohup.out +.vscode/settings.json # Byte-compiled / optimized / DLL files __pycache__/ 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() diff --git a/generate_designs.ipynb b/generate_designs.ipynb deleted file mode 100644 index 4a352e7..0000000 --- a/generate_designs.ipynb +++ /dev/null @@ -1,143 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "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", - "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", - "save_design = False\n", - "plot_design = True\n", - "n_designs = 3\n", - "\n", - "# Define constants\n", - "AEDT_VERSION = \"2024.2\"\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", - " 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", - "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)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "w_mins_base = np.array([3, 2.5, 2.5, 2]) - 1.0\n", - "\n", - "# Loop over the number of designs\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", - "\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", - "\n", - " # Delete the rotor\n", - " design.delete_rotor()\n", - "\n", - " # Potentially plot the desing\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)" - ] - }, - { - "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 -} diff --git a/geometry.py b/geometry.py deleted file mode 100644 index dfd4d78..0000000 --- a/geometry.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np -from scipy.interpolate import CubicSpline - -def rotation_matrix(alpha, rad=False): - if not rad: - 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 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))) - - 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 diff --git a/machine_design/__init__.py b/machine_design/__init__.py new file mode 100644 index 0000000..abd45fb --- /dev/null +++ b/machine_design/__init__.py @@ -0,0 +1,9 @@ +from .design import Design +from .generators import ( + FourStupid, + HacklGenerator_OneLambda, + HacklGenerator_OneLambdaTheta, + HacklGenerator_TwoLambdas, + save_params, +) +from .geometry import analyze_results, plot_barriers, rotate diff --git a/design.py b/machine_design/design.py similarity index 72% rename from design.py rename to machine_design/design.py index d5aee45..8234366 100644 --- a/design.py +++ b/machine_design/design.py @@ -1,15 +1,33 @@ -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=None): + def __init__(self, m2d: Maxwell2d) -> None: self.set_parameters() - if m2d is not None: - self.m2d = m2d + self.m2d = m2d + + @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) + obj = cls(m2d) + obj.create_stator() + obj.save_project(file_name) + return obj + + @classmethod + def load(cls, file_name: str, **kwargs) -> "Design": + desktop = Desktop(**kwargs) + desktop.load_project(file_name) + + m2d = Maxwell2d() + m2d.set_active_design("Design01") + + return cls(m2d) - def set_parameters(self): - #materials + def set_parameters(self) -> None: + # materials self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" # main definitions @@ -20,7 +38,7 @@ def set_parameters(self): "SlotNumber": "36", "SlotType": "3", "DiaShaft": "25mm", - "StackLength": "85mm" + "StackLength": "85mm", } # stator slot self.slot_params = { @@ -40,51 +58,52 @@ def set_parameters(self): "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 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): + 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" @@ -124,7 +143,7 @@ def create_stator(self): name="Shaft", ) - #motion setup + # motion setup m2d.assign_rotate_motion( assignment="Band", coordinate_system="Global", @@ -182,7 +201,7 @@ def create_stator(self): 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], @@ -192,9 +211,7 @@ def create_stator(self): 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 @@ -214,9 +231,7 @@ def create_stator(self): # 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))", @@ -243,16 +258,14 @@ def create_stator(self): ], 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", @@ -279,9 +292,7 @@ def create_stator(self): 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", @@ -308,9 +319,7 @@ def create_stator(self): 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", @@ -337,9 +346,7 @@ def create_stator(self): 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( @@ -357,18 +364,16 @@ def create_stator(self): 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))" @@ -379,7 +384,7 @@ def create_stator(self): setup.update() m2d.validate_simple() - #ooutput variables + # ooutput variables output_vars = { "pos": "(Moving1.Position -InitPos) * Poles/2", "cos0": "cos(pos)", @@ -407,22 +412,25 @@ def create_stator(self): 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", @@ -437,49 +445,43 @@ def create_stator(self): plot_name=v, ) - def add_rotor(self): - rotor_id = self.m2d.modeler.create_polyline( - points=self.rot_points, segment_type=[ "Arc","Line", "Arc"], cover_surface=True, name="Rotor" + 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 rotor_id.material_name = self.Fe 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_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)))) - - 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) + 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" + ) + + # Remove the barrier + 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 + assert m2d.mesh is not None + assert m2d.post is not None m2d.mesh.assign_length_mesh( assignment=self.rotor_id, @@ -488,52 +490,38 @@ def compute(self, NUM_CORES=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" - ) - return solutions.data_magnitude() - - def delete_rotor(self): + 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, **kwargs): - show = kwargs.pop('show', False) - view = kwargs.pop('view', 'xy') + 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): - self.m2d.close_desktop() - def analyze_results(self, Tor): - 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 + def close_project(self) -> None: + self.m2d.close_desktop() - 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') + if not val.endswith("mm"): + raise Exception("val must end with mm") return float(val[:-2]) - - def print_results(self, TorAvg, TorRmsAC, TorRippleRms): - 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/machine_design/generators.py b/machine_design/generators.py new file mode 100644 index 0000000..fd271bd --- /dev/null +++ b/machine_design/generators.py @@ -0,0 +1,362 @@ +import pickle +from abc import ABC, abstractmethod +from typing import Any + +import numpy as np +from design import Design +from geometry import rotate +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: + angles = np.linspace(np.radians(start_deg), np.radians(end_deg), n_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 + + +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: + 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: + 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 + + @abstractmethod + def random_parameters(self) -> Any: + pass + + @abstractmethod + 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__ + + 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] + + # 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) + + 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]): + f = signed_distance(barrier, a, b, c) + f_is_positive = f > 0 + 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]: + 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: + 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__(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] + y1 = [y_min, y_mid, y_max1] + 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) + + 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) + # 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: + 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 + + @property + def bounds(self) -> tuple[np.ndarray, np.ndarray]: + 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 = 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 + 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 X_to_params(self, X: np.ndarray): + return X[:4], X[4], X[5], X[6], X[7] + + 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) + 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.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.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) + + @abstractmethod + 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) + 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 = [] + 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:])) + barriers.append(barrier) + return barriers + + 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) + 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]) + r2 = np.array([y_bezier, x_bezier]) + r3 = np.array([y_end, x_end]) + + # 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 + + 3 * (1 - z_vals) * z_vals**2 * r2 + + 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) + + @property + 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] + 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 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 + + +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) + + @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] + 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 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: + xx, yy = rotate(x, y, self.theta) + return xx[0], yy[0] + return x, y + + +class HacklGenerator_TwoLambdas(AbstractHacklGenerator): + def __init__(self, design, r_stator_end, **kwargs): + 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 + 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] + 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 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]: + lam = self.lam_inner if is_inner else self.lam_outer + return lam * x + (1 - lam) * y, y diff --git a/machine_design/geometry.py b/machine_design/geometry.py new file mode 100644 index 0000000..ecac422 --- /dev/null +++ b/machine_design/geometry.py @@ -0,0 +1,36 @@ +import matplotlib.pyplot as plt +import numpy as np + + +def rotation_matrix(alpha, rad=False): + if not rad: + 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) + 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]) + plt.axis("equal") + if title is not None: + plt.title(title) + if file_name is not 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 + return TorAvg, TorRmsAC, TorRippleRms 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 +} diff --git a/notebooks/run.py b/notebooks/run.py new file mode 100644 index 0000000..b19e2b1 --- /dev/null +++ b/notebooks/run.py @@ -0,0 +1,101 @@ +# export ANSYSEM_ROOT241=/data/AnsysEM/v241/Linux64 + +import os + +import numpy as np +import pandas as pd + +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" +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 = 20 + +# 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: + params = generator.random_parameters() + 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() + + # 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") + save_params(params, f"{file_name}.pkl") + + 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() diff --git a/notebooks/run_optimization.py b/notebooks/run_optimization.py new file mode 100644 index 0000000..d7e28b2 --- /dev/null +++ b/notebooks/run_optimization.py @@ -0,0 +1,202 @@ +# export ANSYSEM_ROOT241=/data/AnsysEM/v241/Linux64 + +import os +import pickle +from collections.abc import Iterable + +import numpy as np +import pandas as pd +import torch +from botorch import fit_gpytorch_mll +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.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 + + +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, + ) + +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_TwoLambdas(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, -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) + 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) + + # Define acquisition function + acq = qLogExpectedHypervolumeImprovement( + model=model, + ref_point=ref_point.tolist(), + partitioning=partitioning, + constraints=[ripple_constraint], + ) + + # 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) + 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() 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 04ea1c0..0000000 --- a/requirements.txt +++ /dev/null @@ -1,4 +0,0 @@ -pyaedt -pyvista[jupyter] -trame -matplotlib \ No newline at end of file