Created
January 4, 2019 03:13
-
-
Save kingoflolz/1c5b6960ab84f0dcc1651a842af39cdf to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import random | |
| import subprocess | |
| import numpy as np | |
| import scipy.optimize | |
| import os | |
| import time | |
| import copy | |
| import multiprocessing | |
| class Simulation: | |
| def __init__(self): | |
| self.commands = ["newdocument(0)", 'mi_probdef(0, "millimeters", "planar", 1E-8)'] | |
| self.commands += ['mi_addboundprop("dirichlet", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)'] | |
| self.points = {} | |
| self.materials = {} | |
| self.circuits = 0 | |
| self.boundaries = 0 | |
| def run_sim(self, id): | |
| self.commands += ['mi_saveas("Untitled{}.fem")'.format(id), "mi_analyse(1)", "mi_loadsolution()"] | |
| def save(self, varname, id): | |
| self.commands += ['f = openfile("femmoutfile{}", "a+")'.format(id), 'write(f, {})'.format(varname), | |
| 'write(f, "\\n")', 'flush()', 'closefile(f)'] | |
| def calculate_force(self, blocks, id): | |
| self.commands += ['mo_hidedensityplot()'] | |
| self.commands += ['mo_hidecontourplot()'] | |
| self.commands += ['mo_seteditmode("area")'] | |
| for i in blocks: | |
| self.commands += ['mo_selectblock({}, {})'.format(i[0], i[1])] | |
| # self.commands += ["pause()"] | |
| self.commands += ['force = mo_blockintegral(18)', | |
| 'mo_close()'] | |
| self.save("force", id) | |
| def zoom_useful(self): | |
| x = [] | |
| y = [] | |
| for k, v in self.points.items(): | |
| x.append(k[0]) | |
| y.append(k[1]) | |
| self.commands += ["mo_zoom({}, {}, {}, {})".format(min(x), min(y), max(x), max(y))] | |
| self.commands += ['mo_showdensityplot(1,0,2,0,"bmag")'] | |
| self.commands += ['mo_showcontourplot(19,-0.01,0.01,real)'] | |
| def take_cap(self, n, id): | |
| self.commands += ['mo_savebitmap("id{}n{}.bmp")'.format(id, n)] | |
| def move_group(self, group, x, y): | |
| self.commands += ["mi_clearselected()"] | |
| # self.commands += ["pause()"] | |
| self.commands += ['mi_selectgroup({})'.format(group)] | |
| # self.commands += ["pause()"] | |
| self.commands += ['mi_movetranslate({}, {}, 4)'.format(x, y)] | |
| # self.commands += ["pause()"] | |
| self.commands += ["mi_clearselected()"] | |
| def add_point(self, x, y, group=0): | |
| if (x, y) not in self.points: | |
| self.commands += ["mi_addnode({}, {})".format(x, y)] | |
| self.points[(x, y)] = 0 | |
| def draw_line(self, x1, y1, x2, y2, group=0, propname=""): | |
| self.add_point(x1, y1) | |
| self.add_point(x2, y2) | |
| self.commands += ["mi_addsegment({}, {}, {}, {})".format(x1, y1, x2, y2)] | |
| self.commands += ["mi_selectsegment({}, {})".format((x1 + x2) / 2, (y1 + y2) / 2)] | |
| self.commands += ['mi_setsegmentprop("{}", 0, 0, 0, {})'.format(propname, group)] | |
| self.commands += ["mi_clearselected()"] | |
| def rect_with_mat(self, x1, y1, x2, y2, name, mag_angle=0, current=0, group=0): | |
| self.draw_contour(0, 0, [ | |
| (x1, y1), | |
| (x1, y2), | |
| (x2, y2), | |
| (x2, y1), | |
| ], group=group) | |
| self.set_materials((x1 + x2) / 2, (y1 + y2) / 2, name, mag_angle, current) | |
| def draw_contour(self, x, y, l, group=0): | |
| for iind, i in enumerate(l): | |
| j = l[iind - 1] | |
| self.draw_line(x + i[0], y + i[1], x + j[0], y + j[1], group=group) | |
| def create_boundaries(self): | |
| x = [] | |
| y = [] | |
| for k, v in self.points.items(): | |
| x.append(k[0]) | |
| y.append(k[1]) | |
| x_center = (max(x) - min(x)) / 2 + min(x) | |
| y_center = (max(y) - min(y)) / 2 + min(y) | |
| r = max(max(x) - min(x), max(y) - min(y)) | |
| self.commands += ["mi_makeABC(7, {}, {}, {}, {})".format(r, x_center, y_center, 0)] | |
| def create_apb(self, x1, y1, x2, y2, x3, y3, x4, y4, group=0): | |
| self.boundaries += 1 | |
| self.commands += ['mi_addboundprop("apb{}", 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0)'.format(self.boundaries)] | |
| self.draw_line(x1, y1, x2, y2, propname="apb{}".format(self.boundaries), group=group) | |
| self.draw_line(x3, y3, x4, y4, propname="apb{}".format(self.boundaries), group=group) | |
| def create_dirichlet_line(self, x1, y1, x2, y2, group=0): | |
| self.draw_line(x1, y1, x2, y2, propname="dirichlet".format(self.boundaries), group=group) | |
| def add_materials(self, name): | |
| if name not in self.materials: | |
| self.materials[name] = 0 | |
| self.commands += ['mi_getmaterial("{}")'.format(name)] | |
| def set_materials(self, x, y, name, mag_angle=0, current=0, group=0): | |
| self.add_materials(name) | |
| current_str = "" | |
| if current != 0: | |
| self.circuits += 1 | |
| self.commands += ['mi_addcircprop("circuit{}", {}, 1)'.format(self.circuits, current, 1)] | |
| current_str = "circuit" + str(self.circuits) | |
| self.commands += ["mi_clearselected()"] | |
| self.commands += ["mi_addblocklabel({}, {})".format(x, y)] | |
| self.commands += ["mi_selectlabel({}, {})".format(x, y)] | |
| self.commands += ['mi_setblockprop("{}", 1, 0, "{}", {}, {}, 1)'.format(name, current_str, mag_angle, group)] | |
| self.commands += ["mi_clearselected()"] | |
| def output(self): | |
| return "\n\r".join(self.commands) | |
| def quit(self): | |
| self.commands += ["quit()"] | |
| class Motor: | |
| def __init__(self, id, stator_dist=9.0): | |
| self.teeth_height = 2 | |
| self.notch_height = 6 | |
| self.teeth_width = 8 | |
| self.notch_width = 2 | |
| self.air_gap = 0.5 | |
| self.back_iron_thickness = 3 | |
| self.stator_dist = stator_dist | |
| self.magnet_width = self.stator_dist * 6 / 7 - 1 | |
| self.magnet_thickness = 2 | |
| self.offset = 0 | |
| self.sim = Simulation() | |
| self.id = id | |
| def draw_outer_magnets(self, x, a): | |
| if a: | |
| a1 = 90 | |
| else: | |
| a1 = 270 | |
| self.sim.rect_with_mat(x, (self.notch_height + self.teeth_height * 2) / 2 + self.air_gap, | |
| x + self.magnet_width, | |
| (self.notch_height + self.teeth_height * 2) / 2 + self.air_gap + self.magnet_thickness, | |
| "NdFeB 40 MGOe", mag_angle=a1) | |
| self.sim.rect_with_mat(x, -(self.notch_height + self.teeth_height * 2) / 2 - self.air_gap, | |
| x + self.magnet_width, | |
| -(self.notch_height + self.teeth_height * 2) / 2 - self.air_gap - self.magnet_thickness, | |
| "NdFeB 40 MGOe", mag_angle=a1) | |
| def draw_stator(self, x, y, cur): | |
| self.sim.draw_contour(x, y, [(self.notch_width / 2, self.notch_height / 2), | |
| (self.notch_width / 2, -self.notch_height / 2), | |
| (self.teeth_width / 2, -self.notch_height / 2), | |
| (self.teeth_width / 2, -self.notch_height / 2 - self.teeth_height), | |
| (-self.teeth_width / 2, -self.notch_height / 2 - self.teeth_height), | |
| (-self.teeth_width / 2, -self.notch_height / 2), | |
| (-self.notch_width / 2, -self.notch_height / 2), | |
| (-self.notch_width / 2, self.notch_height / 2), | |
| (-self.teeth_width / 2, self.notch_height / 2), | |
| (-self.teeth_width / 2, self.notch_height / 2 + self.teeth_height), | |
| (self.teeth_width / 2, self.notch_height / 2 + self.teeth_height), | |
| (self.teeth_width / 2, self.notch_height / 2), | |
| ], group=1) | |
| self.sim.draw_line(self.teeth_width / 2 + x, self.notch_height / 2 + y, self.teeth_width / 2 + x, | |
| -self.notch_height / 2 + y, | |
| group=1) | |
| self.sim.draw_line(-self.teeth_width / 2 + x, self.notch_height / 2 + y, -self.teeth_width / 2 + x, | |
| -self.notch_height / 2 + y, | |
| group=1) | |
| copper_width = (self.teeth_width - self.notch_width) / 2 | |
| current = 50 * (copper_width * self.notch_height) * cur | |
| self.sim.set_materials(x, y, "M-22 Steel", group=1) | |
| self.sim.set_materials(self.notch_width / 2 + 0.01 + x, y, "Copper", current=int(current), group=1) | |
| self.sim.set_materials(-self.notch_width / 2 - 0.01 + x, y, "Copper", current=-int(current), group=1) | |
| def build(self, steps, pic): | |
| for iind, i in enumerate([1, -1, 0.5, -0.5, 0.5, -0.5]): | |
| self.draw_stator(((iind - 2.5) * self.stator_dist), 0, i) | |
| magnet_space = self.stator_dist * 6 / 7 - self.magnet_width | |
| for i in range(7): | |
| self.draw_outer_magnets((i - 3.5) * (self.magnet_width + magnet_space) + magnet_space / 2, i % 2) | |
| min_x = -3.5 * (self.magnet_width + magnet_space) | |
| max_x = 3.5 * (self.magnet_width + magnet_space) | |
| stator_end_y = (self.notch_height + self.teeth_height * 2) / 2 | |
| back_iron_start_y = stator_end_y + self.air_gap + self.magnet_thickness | |
| back_iron_end_y = stator_end_y + self.air_gap + self.magnet_thickness + self.back_iron_thickness | |
| air_end_y = back_iron_end_y + 2 | |
| magnet_start_y = stator_end_y + self.air_gap | |
| self.sim.rect_with_mat(min_x, | |
| back_iron_start_y, | |
| max_x, | |
| back_iron_end_y, | |
| "1018 Steel") | |
| self.sim.rect_with_mat(min_x, | |
| -back_iron_start_y, | |
| max_x, | |
| -back_iron_end_y, | |
| "1018 Steel") | |
| self.sim.set_materials(0, (self.notch_height + self.air_gap / 2) / 2 + self.teeth_height, "Air") | |
| self.sim.set_materials(0, back_iron_end_y + 0.1, "Air") | |
| self.sim.set_materials(0, -back_iron_end_y - 0.1, "Air") | |
| self.sim.create_apb(min_x, back_iron_start_y, min_x, back_iron_end_y, max_x, back_iron_start_y, max_x, | |
| back_iron_end_y) | |
| self.sim.create_apb(min_x, -back_iron_start_y, min_x, -back_iron_end_y, max_x, -back_iron_start_y, max_x, | |
| -back_iron_end_y) | |
| self.sim.create_apb(min_x, magnet_start_y, min_x, stator_end_y, max_x, magnet_start_y, max_x, | |
| stator_end_y) | |
| self.sim.create_apb(min_x, -magnet_start_y, min_x, -stator_end_y, max_x, -magnet_start_y, max_x, | |
| -stator_end_y) | |
| self.sim.create_apb(min_x, magnet_start_y, min_x, back_iron_start_y, max_x, magnet_start_y, max_x, | |
| back_iron_start_y) | |
| self.sim.create_apb(min_x, -magnet_start_y, min_x, -back_iron_start_y, max_x, -magnet_start_y, max_x, | |
| -back_iron_start_y) | |
| self.sim.create_apb(min_x, air_end_y, min_x, back_iron_end_y, max_x, air_end_y, max_x, back_iron_end_y) | |
| self.sim.create_apb(min_x, -air_end_y, min_x, -back_iron_end_y, max_x, -air_end_y, max_x, -back_iron_end_y) | |
| self.sim.create_dirichlet_line(min_x, air_end_y, max_x, air_end_y) | |
| self.sim.create_dirichlet_line(min_x, -air_end_y, max_x, -air_end_y) | |
| self.sim.create_apb(min_x, stator_end_y, min_x, -stator_end_y, max_x, stator_end_y, max_x, -stator_end_y, | |
| group=1) | |
| b = 0 | |
| m = (magnet_space + self.magnet_width) / steps | |
| self.sim.move_group(1, b, 0) | |
| for i in range(steps): | |
| k = b + m * i | |
| self.sim.run_sim(self.id) | |
| if pic: | |
| print("taking image") | |
| self.sim.zoom_useful() | |
| self.sim.take_cap(i, self.id) | |
| force_blocks = [] | |
| for i in range(6): | |
| o = (i - 2.5) * self.stator_dist | |
| force_blocks += [(o + k, 0), (o + k + self.teeth_width / 2 - 0.01, 0), | |
| (o + k - self.teeth_width / 2 + 0.01, 0)] | |
| self.sim.calculate_force(force_blocks, self.id) | |
| self.sim.move_group(1, m, 0) | |
| self.sim.quit() | |
| open("sim{}.lua".format(self.id), mode="w+").write(self.sim.output()) | |
| def run(self): | |
| a = subprocess.Popen(["/home/ben/.wine/drive_c/femm42/bin/femm.exe", | |
| "-lua-script=/home/ben/PycharmProjects/axial_motor2/sim{}.lua".format(self.id)]) | |
| return a | |
| def analyze(self): | |
| torques = [] | |
| for line in open("femmoutfile{}".format(self.id), mode="r+").readlines(): | |
| torques.append(abs(float(line.strip()))) | |
| os.remove("femmoutfile{}".format(self.id)) | |
| os.remove("Untitled{}.fem".format(self.id)) | |
| os.remove("sim{}.lua".format(self.id)) | |
| os.remove("Untitled{}.ans".format(self.id)) | |
| return torques | |
| def ef(f, epsilon, xk, f0, k, args): | |
| ei = np.zeros((len(xk),), float) | |
| ei[k] = 1.0 | |
| d = epsilon * ei | |
| return (f(*((xk + d,) + args)) - f0) / d[k] | |
| if __name__ == '__main__': | |
| def fitness(vars, steps=1, pic=False): | |
| start = time.time() | |
| r = int(time.time() * 1000) * 1000 + random.randint(0, 1000) | |
| time.sleep(0.01) | |
| m1 = Motor(r) | |
| m1.teeth_height = vars[0] | |
| m1.notch_height = vars[1] | |
| m1.notch_width = vars[2] | |
| m1.air_gap = vars[3] | |
| m1.back_iron_thickness = vars[4] | |
| m1.magnet_thickness = vars[5] | |
| m2 = copy.deepcopy(m1) | |
| m2.id = int(time.time() * 1000) * 1000 + random.randint(0, 1000) | |
| time.sleep(0.01) | |
| m2.stator_dist = 12.251 | |
| m3 = copy.deepcopy(m2) | |
| m3.id = int(time.time() * 1000) * 1000 + random.randint(0, 1000) | |
| time.sleep(0.01) | |
| m3.teeth_width = 11.215 | |
| m3.notch_width += 11.215 - 8 | |
| m3.magnet_width = 6 * 11.215 / 7 - 1 | |
| m4 = copy.deepcopy(m3) | |
| m4.stator_dist = 15.24 | |
| m4.id = int(time.time() * 1000) * 1000 + random.randint(0, 1000) | |
| time.sleep(0.01) | |
| m4.stator_dist = 15.248 | |
| m1.build(steps, pic) | |
| m2.build(steps, pic) | |
| m3.build(steps, pic) | |
| m4.build(steps, pic) | |
| a = m1.run() | |
| b = m2.run() | |
| c = m3.run() | |
| d = m4.run() | |
| a.wait() | |
| b.wait() | |
| c.wait() | |
| d.wait() | |
| t1 = np.array(m1.analyze()) | |
| t2 = np.array(m2.analyze()) | |
| t3 = np.array(m3.analyze()) | |
| t4 = np.array(m4.analyze()) | |
| print(t1) | |
| print(t2) | |
| print(t3) | |
| print(t4) | |
| torque = 2 * (t1 * 3 * 0.0185 + t2 * 3 * 0.0215 + t3 * 2.8 * 0.0245 + t4 * 2.8 * 0.0245) | |
| actual = -torque.max() / ( | |
| 2 * m1.teeth_height + m1.notch_height + 2 * m1.air_gap + 2 * m1.back_iron_thickness + 2 * m1.magnet_thickness + 20) | |
| print(vars) | |
| print(actual) | |
| print(torque) | |
| print("thickness: {}".format( | |
| 2 * m1.teeth_height + m1.notch_height + 2 * m1.air_gap + 2 * m1.back_iron_thickness + 2 * m1.magnet_thickness + 20)) | |
| print("simulation took {} s".format(time.time() - start)) | |
| return actual | |
| bounds = [(0, 8), (1, 15), (0.5, 3.5), (0.5, 2), (0.5, 3), (1, 5)] | |
| constraints = [ | |
| scipy.optimize.LinearConstraint(np.array([0, 0, 1, 0, 0, 0]), 0.1, 4) | |
| ] | |
| pool = multiprocessing.Pool(2) | |
| def jac(xk, epsilon=0.01, args=(), f0=None): | |
| """ | |
| See ``approx_fprime``. An optional initial function value arg is added. | |
| """ | |
| if f0 is None: | |
| f0 = fitness(*((xk,) + args)) | |
| grad = pool.starmap(ef, [(fitness, epsilon, xk, f0, i, args) for i in range(len(xk))]) | |
| print(grad) | |
| return grad | |
| start = np.array([0.8092011 , 8.52388944, 2.39536697, 0.5, 2.43616628,2.37874529]) | |
| fitness(start, steps=1, pic=True) | |
| print(scipy.optimize.minimize(fitness, start, bounds=bounds, | |
| constraints=constraints, | |
| options={'eps': 0.01, "maxiter": 200}, jac=jac)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment