geneticOptimizationSliderCrank.py
You can view and download this file on Github: geneticOptimizationSliderCrank.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Slider crank model with verification in MATLAB for machine dynamics course
5# optionally, the slider crank is mounted on a floating frame, leading to vibrations
6# if the system is unbalanced
7# Use this example in combination with cmd: 'python resultsMonitor.py solution/geneticSliderCrank.txt'
8#
9# Author: Johannes Gerstmayr
10# Date: 2019-12-07 (created)
11# 2021-01-10 (adapted for genetic optimization)
12#
13# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
14#
15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16
17import exudyn as exu
18from exudyn.itemInterface import *
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import exudyn.graphics as graphics #only import if it does not conflict
21from exudyn.processing import GeneticOptimization, ParameterVariation, PlotOptimizationResults2D
22
23import numpy as np #for postprocessing
24import os
25from time import sleep
26
27useGraphics = False
28L1=0.1
29L2=0.3
30m1=0.4
31m2=0.2
32m3=0.1
33s1opt = -L1*(m2+m3)/m1 #-0.075
34s2opt = -m3/m2*L2 #-0.15
35
36#this is the function which is repeatedly called from ParameterVariation
37#parameterSet contains dictinary with varied parameters
38def ParameterFunction(parameterSet):
39 SC = exu.SystemContainer()
40 mbs = SC.AddSystem()
41
42 #++++++++++++++++++++++++++++++++++++++++++++++
43 #++++++++++++++++++++++++++++++++++++++++++++++
44 #store default parameters in structure (all these parameters can be varied!)
45 class P: pass #create emtpy structure for parameters; simplifies way to update parameters
46 P.s1=L1*0.5
47 P.s2=L2*0.5
48 P.h=0.002
49 P.computationIndex = ''
50
51 # #now update parameters with parameterSet (will work with any parameters in structure P)
52 for key,value in parameterSet.items():
53 setattr(P,key,value)
54
55 globalSettings = parameterSet['functionData']
56 stiffness = globalSettings['stiffness']
57 computationIndex = parameterSet['computationIndex'] #could be used for temporary variables when using multithreading
58
59 #++++++++++++++++++++++++++++++++++++++++++++++
60 #++++++++++++++++++++++++++++++++++++++++++++++
61 #START HERE: create parameterized model
62
63 testCases = 1 #floating body
64 nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
65 mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
66
67
68 #++++++++++++++++++++++++++++++++
69 #floating body to mount slider-crank mechanism
70 constrainGroundBody = (testCases == 0) #use this flag to fix ground body
71
72 #graphics for floating frame:
73 gFloating = graphics.BrickXYZ(-0.25, -0.25, -0.1, 0.8, 0.25, -0.05, color=[0.3,0.3,0.3,1.])
74
75 if constrainGroundBody:
76 floatingRB = mbs.AddObject(ObjectGround(referencePosition=[0,0,0], visualization=VObjectGround(graphicsData=[gFloating])))
77 mFloatingN = mbs.AddMarker(MarkerBodyPosition(bodyNumber = floatingRB, localPosition=[0,0,0]))
78 else:
79 nFloating = mbs.AddNode(Rigid2D(referenceCoordinates=[0,0,0], initialVelocities=[0,0,0]));
80 mFloatingN = mbs.AddMarker(MarkerNodePosition(nodeNumber=nFloating))
81 floatingRB = mbs.AddObject(RigidBody2D(physicsMass=2, physicsInertia=1, nodeNumber=nFloating, visualization=VObjectRigidBody2D(graphicsData=[gFloating])))
82 mRB0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=0))
83 mRB1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=1))
84 mRB2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=2))
85
86 #add spring dampers for reference frame:
87 k=stiffness #stiffness of floating body
88 d=k*0.01
89 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB0], stiffness=k, damping=d))
90 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB1], stiffness=k, damping=d))
91 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB2], stiffness=k, damping=d))
92 mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mRB2]))
93
94
95
96 #++++++++++++++++++++++++++++++++
97 #nodes and bodies
98 omega=2*pi/60*300 #3000 rpm
99 M=0.1 #torque (default: 0.1)
100
101 s1L=-P.s1
102 s1R=L1-P.s1
103 s2L=-P.s2
104 s2R=L2-P.s2
105
106 #lambda=L1/L2
107 J1=(m1/12.)*L1**2 #inertia w.r.t. center of mass
108 J2=(m2/12.)*L2**2 #inertia w.r.t. center of mass
109
110 ty = 0.05 #thickness
111 tz = 0.05 #thickness
112
113 graphics1 = graphics.RigidLink(p0=[s1L,0,-0.5*tz],p1=[s1R,0,-0.5*tz],
114 axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
115 thickness=0.8*ty, width=[tz,tz], color=graphics.color.steelblue,nTiles=16)
116
117 graphics2 = graphics.RigidLink(p0=[s2L,0,0.5*tz],p1=[s2R,0,0.5*tz],
118 axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
119 thickness=0.8*ty, width=[tz,tz], color=graphics.color.lightred,nTiles=16)
120
121 #crank:
122 nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[P.s1,0,0],
123 initialVelocities=[0,0,0]));
124 oRigid1 = mbs.AddObject(RigidBody2D(physicsMass=m1,
125 physicsInertia=J1,
126 nodeNumber=nRigid1,
127 visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
128
129 #connecting rod:
130 nRigid2 = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+P.s2,0,0],
131 initialVelocities=[0,0,0]));
132 oRigid2 = mbs.AddObject(RigidBody2D(physicsMass=m2,
133 physicsInertia=J2,
134 nodeNumber=nRigid2,
135 visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
136
137
138 #++++++++++++++++++++++++++++++++
139 #slider:
140 c=0.025 #dimension of mass
141 graphics3 = graphics.BrickXYZ(-c,-c,-c*2,c,c,0,graphics.color.grey)
142
143 #nMass = mbs.AddNode(Point2D(referenceCoordinates=[L1+L2,0]))
144 #oMass = mbs.AddObject(MassPoint2D(physicsMass=m3, nodeNumber=nMass,visualization=VObjectMassPoint2D(graphicsData= [graphics3])))
145 nMass = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+L2,0,0]))
146 oMass = mbs.AddObject(RigidBody2D(physicsMass=m3, physicsInertia=0.001*m3, nodeNumber=nMass,visualization=VObjectRigidBody2D(graphicsData= [graphics3])))
147
148 #++++++++++++++++++++++++++++++++
149 #markers for joints:
150 mR1Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRigid1, localPosition= [s1L,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
151 mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[s1R,0.,0.])) #end point; connection to connecting rod
152
153 mR2Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition= [s2L,0.,0.])) #connection to crank
154 mR2Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition=[s2R,0.,0.])) #end point; connection to slider
155
156 mMass = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[ 0.,0.,0.]))
157 mG0 = mFloatingN
158
159 #++++++++++++++++++++++++++++++++
160 #joints:
161 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1Left]))
162 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mR2Left]))
163 mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR2Right,mMass]))
164
165
166 #prismatic joint:
167 mRigidGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber = floatingRB, localPosition = [L1+L2,0,0]))
168 mRigidSlider = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oMass, localPosition = [0,0,0]))
169
170 mbs.AddObject(PrismaticJoint2D(markerNumbers=[mRigidGround,mRigidSlider], constrainRotation=True))
171
172
173 #user function for load; switch off load after 1 second
174 userLoadOn = True
175 def userLoad(mbs, t, load):
176 setLoad = 0
177 if userLoadOn:
178 setLoad = load
179 omega = mbs.GetNodeOutput(nRigid1,variableType = exu.OutputVariableType.AngularVelocity)[2]
180 if omega > 2*pi*2:
181 #print("t=",t)
182 userLoadOn = False
183 return setLoad
184
185 #loads and driving forces:
186 mRigid1CoordinateTheta = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid1, coordinate=2)) #angle coordinate is constrained
187 #mbs.AddLoad(LoadCoordinate(markerNumber=mRigid1CoordinateTheta, load = M, loadUserFunction=userLoad)) #torque at crank
188 mbs.AddLoad(LoadCoordinate(markerNumber=mRigid1CoordinateTheta, load = M)) #torque at crank
189
190 #write motion of support frame:
191 sFloating = mbs.AddSensor(SensorNode(nodeNumber=nFloating,
192 storeInternal=True,
193 outputVariableType=exu.OutputVariableType.Position))
194
195 #++++++++++++++++++++++++++++++++
196 #assemble, adjust settings and start time integration
197 mbs.Assemble()
198
199 simulationSettings = exu.SimulationSettings() #takes currently set values or default values
200 tEnd = 3
201
202 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/P.h)
203 simulationSettings.timeIntegration.endTime = tEnd
204
205
206 simulationSettings.solutionSettings.solutionWritePeriod = 2e-3
207 simulationSettings.solutionSettings.writeSolutionToFile = useGraphics
208
209 simulationSettings.timeIntegration.newton.useModifiedNewton = True
210 simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
211 simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
212
213 #++++++++++++++++++++++++++++++++++++++++++
214 #solve index 2 / trapezoidal rule:
215 simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
216 simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
217
218 dSize = 0.02
219 SC.visualizationSettings.nodes.defaultSize = dSize
220 SC.visualizationSettings.markers.defaultSize = dSize
221 SC.visualizationSettings.bodies.defaultSize = [dSize, dSize, dSize]
222 SC.visualizationSettings.connectors.defaultSize = dSize
223
224 #data obtained from SC.renderer.GetState(); use np.round(d['modelRotation'],4)
225 SC.visualizationSettings.openGL.initialModelRotation = [[ 0.87758, 0.04786, -0.47703],
226 [ 0. , 0.995 , 0.09983],
227 [ 0.47943, -0.08761, 0.8732]]
228 SC.visualizationSettings.openGL.initialZoom = 0.47
229 SC.visualizationSettings.openGL.initialCenterPoint = [0.192, -0.0039,-0.075]
230 SC.visualizationSettings.openGL.initialMaxSceneSize = 0.4
231 SC.visualizationSettings.general.autoFitScene = False
232 #SC.renderer.DoIdleTasks()
233
234 if useGraphics:
235 SC.renderer.Start()
236
237 mbs.SolveDynamic(simulationSettings)
238
239 if useGraphics:
240 SC.renderer.DoIdleTasks()
241 SC.renderer.Stop() #safely close rendering window!
242
243 #++++++++++++++++++++++++++++++++++++++++++
244 #evaluate error:
245 #data = np.loadtxt(sensorFileName, comments='#', delimiter=',')
246 data = mbs.GetSensorStoredData(sFloating)
247
248 errorNorm = max(abs(data[:,1])) + max(abs(data[:,2])) #max displacement in x and y direction
249
250 if useGraphics:
251 print("max. oszillation=", errorNorm)
252
253 mbs.PlotSensor(sensorNumbers=[sFloating,sFloating], components=[0,1])
254
255 del mbs
256 del SC
257
258 return errorNorm
259 #++++++++++++++++++++++++++++++++++++++++++
260
261import matplotlib.pyplot as plt
262import matplotlib.ticker as ticker
263
264doOptimize = True
265#now perform parameter variation
266if __name__ == '__main__': #include this to enable parallel processing
267 if doOptimize:
268 import time
269
270 #some data which shall not be optimized, but passed to objectiveFunction (e.g. in variation calculations)
271 functionData = {'stiffness':5000}
272
273 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
274 #GeneticOptimization
275 start_time = time.time()
276 [pOpt, vOpt, pList, values] = GeneticOptimization(objectiveFunction = ParameterFunction,
277 parameters = {'s1':(-L1,L1), 's2':(-L2,L2)}, #parameters provide search range
278 parameterFunctionData = functionData,
279 numberOfGenerations = 30,
280 populationSize = 50,
281 elitistRatio = 0.1,
282 crossoverProbability = 0.1,
283 rangeReductionFactor = 0.5,
284 addComputationIndex=True,
285 randomizerInitialization=0, #for reproducible results
286 #distanceFactor = 0.1, #for this example only one significant minimum
287 debugMode=False,
288 useMultiProcessing=True, #may be problematic for test
289 showProgress=True,
290 resultsFile = 'solution/geneticSliderCrank.txt',
291 )
292 #exu.Print("--- %s seconds ---" % (time.time() - start_time))
293
294 exu.Print("[pOpt, vOpt]=", [pOpt, vOpt])
295 u = vOpt
296 exu.Print("optimum=",u)
297 # using files:
298 # [pOpt, vOpt]= [{'s1': -0.07497827333782427, 's2': -0.14943029494085874}, 3.4312580948e-05]
299 # optimum= 3.4312580948e-05
300
301 # using internal storage:
302 # [pOpt, vOpt]= [{'s1': -0.07497827333782427, 's2': -0.14943029494085874}, 3.431258094752888e-05]
303 # optimum= 3.431258094752888e-05
304
305 if False:
306 # from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
307 import matplotlib.pyplot as plt
308
309 plt.close('all')
310 [figList, axList] = PlotOptimizationResults2D(pList, values, yLogScale=True)
311 else:
312 useGraphics = True
313 parameterSet = {'s1':L1*0.5, 's2':L2*0.5, 'h':1e-5}
314 #parameterSet = {'s1':-0.075, 's2':-0.15, 'h':1e-5}
315 ParameterFunction(parameterSet)