You can view and download this file on Github: geneticOptimizationSliderCrank.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Slider crank model with verification in MATLAB for machine dynamics course
# optionally, the slider crank is mounted on a floating frame, leading to vibrations
# if the system is unbalanced
# Use this example in combination with cmd: 'python resultsMonitor.py solution/geneticSliderCrank.txt'
#
# Author: Johannes Gerstmayr
# Date: 2019-12-07 (created)
# 2021-01-10 (adapted for genetic optimization)
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.itemInterface import *
from exudyn.utilities import *
from exudyn.processing import GeneticOptimization, ParameterVariation, PlotOptimizationResults2D
import numpy as np #for postprocessing
import os
from time import sleep
useGraphics = False
L1=0.1
L2=0.3
m1=0.4
m2=0.2
m3=0.1
s1opt = -L1*(m2+m3)/m1 #-0.075
s2opt = -m3/m2*L2 #-0.15
#this is the function which is repeatedly called from ParameterVariation
#parameterSet contains dictinary with varied parameters
def ParameterFunction(parameterSet):
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#++++++++++++++++++++++++++++++++++++++++++++++
#++++++++++++++++++++++++++++++++++++++++++++++
#store default parameters in structure (all these parameters can be varied!)
class P: pass #create emtpy structure for parameters; simplifies way to update parameters
P.s1=L1*0.5
P.s2=L2*0.5
P.h=0.002
P.computationIndex = ''
# #now update parameters with parameterSet (will work with any parameters in structure P)
for key,value in parameterSet.items():
setattr(P,key,value)
#++++++++++++++++++++++++++++++++++++++++++++++
#++++++++++++++++++++++++++++++++++++++++++++++
#START HERE: create parameterized model
testCases = 1 #floating body
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
#++++++++++++++++++++++++++++++++
#floating body to mount slider-crank mechanism
constrainGroundBody = (testCases == 0) #use this flag to fix ground body
#graphics for floating frame:
gFloating = GraphicsDataOrthoCube(-0.25, -0.25, -0.1, 0.8, 0.25, -0.05, color=[0.3,0.3,0.3,1.])
if constrainGroundBody:
floatingRB = mbs.AddObject(ObjectGround(referencePosition=[0,0,0], visualization=VObjectGround(graphicsData=[gFloating])))
mFloatingN = mbs.AddMarker(MarkerBodyPosition(bodyNumber = floatingRB, localPosition=[0,0,0]))
else:
nFloating = mbs.AddNode(Rigid2D(referenceCoordinates=[0,0,0], initialVelocities=[0,0,0]));
mFloatingN = mbs.AddMarker(MarkerNodePosition(nodeNumber=nFloating))
floatingRB = mbs.AddObject(RigidBody2D(physicsMass=2, physicsInertia=1, nodeNumber=nFloating, visualization=VObjectRigidBody2D(graphicsData=[gFloating])))
mRB0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=0))
mRB1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=1))
mRB2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nFloating, coordinate=2))
#add spring dampers for reference frame:
k=5000 #stiffness of floating body
d=k*0.01
mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB0], stiffness=k, damping=d))
mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB1], stiffness=k, damping=d))
mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGround,mRB2], stiffness=k, damping=d))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mRB2]))
#++++++++++++++++++++++++++++++++
#nodes and bodies
omega=2*pi/60*300 #3000 rpm
M=0.1 #torque (default: 0.1)
s1L=-P.s1
s1R=L1-P.s1
s2L=-P.s2
s2R=L2-P.s2
#lambda=L1/L2
J1=(m1/12.)*L1**2 #inertia w.r.t. center of mass
J2=(m2/12.)*L2**2 #inertia w.r.t. center of mass
ty = 0.05 #thickness
tz = 0.05 #thickness
graphics1 = GraphicsDataRigidLink(p0=[s1L,0,-0.5*tz],p1=[s1R,0,-0.5*tz],
axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
thickness=0.8*ty, width=[tz,tz], color=color4steelblue,nTiles=16)
graphics2 = GraphicsDataRigidLink(p0=[s2L,0,0.5*tz],p1=[s2R,0,0.5*tz],
axis0=[0,0,1], axis1=[0,0,1],radius=[0.5*ty,0.5*ty],
thickness=0.8*ty, width=[tz,tz], color=color4lightred,nTiles=16)
#crank:
nRigid1 = mbs.AddNode(Rigid2D(referenceCoordinates=[P.s1,0,0],
initialVelocities=[0,0,0]));
oRigid1 = mbs.AddObject(RigidBody2D(physicsMass=m1,
physicsInertia=J1,
nodeNumber=nRigid1,
visualization=VObjectRigidBody2D(graphicsData= [graphics1])))
#connecting rod:
nRigid2 = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+P.s2,0,0],
initialVelocities=[0,0,0]));
oRigid2 = mbs.AddObject(RigidBody2D(physicsMass=m2,
physicsInertia=J2,
nodeNumber=nRigid2,
visualization=VObjectRigidBody2D(graphicsData= [graphics2])))
#++++++++++++++++++++++++++++++++
#slider:
c=0.025 #dimension of mass
graphics3 = GraphicsDataOrthoCube(-c,-c,-c*2,c,c,0,color4grey)
#nMass = mbs.AddNode(Point2D(referenceCoordinates=[L1+L2,0]))
#oMass = mbs.AddObject(MassPoint2D(physicsMass=m3, nodeNumber=nMass,visualization=VObjectMassPoint2D(graphicsData= [graphics3])))
nMass = mbs.AddNode(Rigid2D(referenceCoordinates=[L1+L2,0,0]))
oMass = mbs.AddObject(RigidBody2D(physicsMass=m3, physicsInertia=0.001*m3, nodeNumber=nMass,visualization=VObjectRigidBody2D(graphicsData= [graphics3])))
#++++++++++++++++++++++++++++++++
#markers for joints:
mR1Left = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oRigid1, localPosition= [s1L,0.,0.])) #support point # MUST be a rigidBodyMarker, because a torque is applied
mR1Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid1, localPosition=[s1R,0.,0.])) #end point; connection to connecting rod
mR2Left = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition= [s2L,0.,0.])) #connection to crank
mR2Right = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid2, localPosition=[s2R,0.,0.])) #end point; connection to slider
mMass = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oMass, localPosition=[ 0.,0.,0.]))
mG0 = mFloatingN
#++++++++++++++++++++++++++++++++
#joints:
mbs.AddObject(RevoluteJoint2D(markerNumbers=[mG0,mR1Left]))
mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR1Right,mR2Left]))
mbs.AddObject(RevoluteJoint2D(markerNumbers=[mR2Right,mMass]))
#prismatic joint:
mRigidGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber = floatingRB, localPosition = [L1+L2,0,0]))
mRigidSlider = mbs.AddMarker(MarkerBodyRigid(bodyNumber = oMass, localPosition = [0,0,0]))
mbs.AddObject(PrismaticJoint2D(markerNumbers=[mRigidGround,mRigidSlider], constrainRotation=True))
#user function for load; switch off load after 1 second
userLoadOn = True
def userLoad(mbs, t, load):
setLoad = 0
if userLoadOn:
setLoad = load
omega = mbs.GetNodeOutput(nRigid1,variableType = exu.OutputVariableType.AngularVelocity)[2]
if omega > 2*pi*2:
#print("t=",t)
userLoadOn = False
return setLoad
#loads and driving forces:
mRigid1CoordinateTheta = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid1, coordinate=2)) #angle coordinate is constrained
#mbs.AddLoad(LoadCoordinate(markerNumber=mRigid1CoordinateTheta, load = M, loadUserFunction=userLoad)) #torque at crank
mbs.AddLoad(LoadCoordinate(markerNumber=mRigid1CoordinateTheta, load = M)) #torque at crank
#write motion of support frame:
sFloating = mbs.AddSensor(SensorNode(nodeNumber=nFloating,
storeInternal=True,
outputVariableType=exu.OutputVariableType.Position))
#++++++++++++++++++++++++++++++++
#assemble, adjust settings and start time integration
mbs.Assemble()
simulationSettings = exu.SimulationSettings() #takes currently set values or default values
tEnd = 3
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/P.h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.solutionSettings.solutionWritePeriod = 2e-3
simulationSettings.solutionSettings.writeSolutionToFile = useGraphics
simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
#++++++++++++++++++++++++++++++++++++++++++
#solve index 2 / trapezoidal rule:
simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
dSize = 0.02
SC.visualizationSettings.nodes.defaultSize = dSize
SC.visualizationSettings.markers.defaultSize = dSize
SC.visualizationSettings.bodies.defaultSize = [dSize, dSize, dSize]
SC.visualizationSettings.connectors.defaultSize = dSize
#data obtained from SC.GetRenderState(); use np.round(d['modelRotation'],4)
SC.visualizationSettings.openGL.initialModelRotation = [[ 0.87758, 0.04786, -0.47703],
[ 0. , 0.995 , 0.09983],
[ 0.47943, -0.08761, 0.8732]]
SC.visualizationSettings.openGL.initialZoom = 0.47
SC.visualizationSettings.openGL.initialCenterPoint = [0.192, -0.0039,-0.075]
SC.visualizationSettings.openGL.initialMaxSceneSize = 0.4
SC.visualizationSettings.general.autoFitScene = False
#mbs.WaitForUserToContinue()
if useGraphics:
exu.StartRenderer()
mbs.SolveDynamic(simulationSettings)
if useGraphics:
SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
#++++++++++++++++++++++++++++++++++++++++++
#evaluate error:
#data = np.loadtxt(sensorFileName, comments='#', delimiter=',')
data = mbs.GetSensorStoredData(sFloating)
errorNorm = max(abs(data[:,1])) + max(abs(data[:,2])) #max displacement in x and y direction
if useGraphics:
print("max. oszillation=", errorNorm)
mbs.PlotSensor(sensorNumbers=[sFloating,sFloating], components=[0,1])
del mbs
del SC
return errorNorm
#++++++++++++++++++++++++++++++++++++++++++
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
doOptimize = True
#now perform parameter variation
if __name__ == '__main__': #include this to enable parallel processing
if doOptimize:
import time
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
#GeneticOptimization
start_time = time.time()
[pOpt, vOpt, pList, values] = GeneticOptimization(objectiveFunction = ParameterFunction,
parameters = {'s1':(-L1,L1), 's2':(-L2,L2)}, #parameters provide search range
numberOfGenerations = 30,
populationSize = 50,
elitistRatio = 0.1,
crossoverProbability = 0.1,
rangeReductionFactor = 0.5,
addComputationIndex=True,
randomizerInitialization=0, #for reproducible results
#distanceFactor = 0.1, #for this example only one significant minimum
debugMode=False,
useMultiProcessing=True, #may be problematic for test
showProgress=True,
resultsFile = 'solution/geneticSliderCrank.txt',
)
#exu.Print("--- %s seconds ---" % (time.time() - start_time))
exu.Print("[pOpt, vOpt]=", [pOpt, vOpt])
u = vOpt
exu.Print("optimum=",u)
# using files:
# [pOpt, vOpt]= [{'s1': -0.07497827333782427, 's2': -0.14943029494085874}, 3.4312580948e-05]
# optimum= 3.4312580948e-05
# using internal storage:
# [pOpt, vOpt]= [{'s1': -0.07497827333782427, 's2': -0.14943029494085874}, 3.431258094752888e-05]
# optimum= 3.431258094752888e-05
if False:
# from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
import matplotlib.pyplot as plt
plt.close('all')
[figList, axList] = PlotOptimizationResults2D(pList, values, yLogScale=True)
else:
useGraphics = True
parameterSet = {'s1':L1*0.5, 's2':L2*0.5, 'h':1e-5}
#parameterSet = {'s1':-0.075, 's2':-0.15, 'h':1e-5}
ParameterFunction(parameterSet)