geneticOptimizationTest.py
You can view and download this file on Github: geneticOptimizationTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: This example performs a genetic algorithm to optimization using a simple
5# mass-spring-damper system; varying mass, spring, ...
6# The objective function is the error compared to
7# a reference solution using reference/nominal values (which are known here, but could originate from a measurement)
8#
9# Author: Johannes Gerstmayr
10# Date: 2020-11-18
11#
12# 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.
13#
14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15
16import exudyn as exu
17from exudyn.itemInterface import *
18from exudyn.processing import GeneticOptimization, ParameterVariation, PlotOptimizationResults2D
19
20import numpy as np #for postprocessing
21import os
22from time import sleep
23
24useGraphics = True #without test
25#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
26#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
27try: #only if called from test suite
28 from modelUnitTests import exudynTestGlobals #for globally storing test results
29 useGraphics = exudynTestGlobals.useGraphics
30except:
31 class ExudynTestGlobals:
32 pass
33 exudynTestGlobals = ExudynTestGlobals()
34#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35
36dataRef = None
37#this is the function which is repeatedly called from ParameterVariation
38#parameterSet contains dictinary with varied parameters
39def ParameterFunction(parameterSet):
40 SC = exu.SystemContainer()
41 mbs = SC.AddSystem()
42 global dataRef #can be avoided, if reference solution is written to file!
43
44
45 #++++++++++++++++++++++++++++++++++++++++++++++
46 #++++++++++++++++++++++++++++++++++++++++++++++
47 #store default parameters in structure (all these parameters can be varied!)
48 class P: pass #create emtpy structure for parameters; simplifies way to update parameters
49
50 P.mass = 1.6 #mass in kg
51 P.spring = 4000 #stiffness of spring-damper in N/m
52 P.damper = 8 #old: 8; damping constant in N/(m/s)
53 P.u0=-0.08 #initial displacement
54 P.v0=1 #initial velocity
55 P.force =80 #force applied to mass
56 P.computationIndex = ''
57
58 #now update parameters with parameterSet (will work with any parameters in structure P)
59 for key,value in parameterSet.items():
60 setattr(P,key,value)
61
62 #++++++++++++++++++++++++++++++++++++++++++++++
63 #++++++++++++++++++++++++++++++++++++++++++++++
64 #create model using parameters P starting here:
65
66 L=0.5 #spring length (for drawing)
67 n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0],
68 initialCoordinates = [P.u0,0,0],
69 initialVelocities= [P.v0,0,0]))
70
71 #ground node
72 nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
73
74 #add mass point (this is a 3D object with 3 coordinates):
75 massPoint = mbs.AddObject(MassPoint(physicsMass = P.mass, nodeNumber = n1))
76
77 #marker for ground (=fixed):
78 groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
79 #marker for springDamper for first (x-)coordinate:
80 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
81
82 #spring-damper between two marker coordinates
83 nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
84 stiffness = P.spring, damping = P.damper))
85
86 #add load:
87 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker, load = P.force))
88 #add sensor:
89
90 sDisp = mbs.AddSensor(SensorObject(objectNumber=nC,
91 storeInternal=True,
92 outputVariableType=exu.OutputVariableType.Displacement))
93
94 #++++++++++++++++++++++++++++++++++++++++++++++
95 #assemble and compute solution
96 mbs.Assemble()
97
98 steps = 100 #number of steps to show solution
99 tEnd = 1 #end time of simulation
100
101 simulationSettings = exu.SimulationSettings()
102 simulationSettings.solutionSettings.writeSolutionToFile = False
103 simulationSettings.solutionSettings.sensorsWritePeriod = 2e-3 #output interval of sensors
104 simulationSettings.timeIntegration.numberOfSteps = steps
105 simulationSettings.timeIntegration.endTime = tEnd
106
107 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping
108
109 mbs.SolveDynamic(simulationSettings)
110
111 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
112 #evaluate difference between reference and optimized solution
113 #reference solution:
114 data = mbs.GetSensorStoredData(sDisp)
115
116 if P.computationIndex=='Ref':
117 dataRef = data
118
119 diff = data[:,1]-dataRef[:,1]
120
121 errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd
122
123 return errorNorm
124
125#generate reference data, also in multi-threaded version this will be computed for all threads!
126refval = ParameterFunction({'computationIndex':'Ref'}) # compute reference solution
127
128#now perform parameter variation
129if __name__ == '__main__': #include this to enable parallel processing
130 import time
131
132 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
133 #GeneticOptimization
134 start_time = time.time()
135 [pOpt, vOpt, pList, values] = GeneticOptimization(objectiveFunction = ParameterFunction,
136 parameters = {'mass':(1,10), 'spring':(100,10000), 'force':(1,1000)}, #parameters provide search range
137 numberOfGenerations = 2,
138 populationSize = 10,
139 elitistRatio = 0.1,
140 crossoverProbability = 0.1,
141 rangeReductionFactor = 0.7,
142 addComputationIndex=True,
143 randomizerInitialization=0, #for reproducible results
144 distanceFactor = 0.1, #for this example only one significant minimum
145 debugMode=False,
146 useMultiProcessing=False, #may be problematic for test
147 showProgress=False,
148 resultsFile = 'solution/geneticOptimizationTest.txt',
149 )
150 exu.Print("--- %s seconds ---" % (time.time() - start_time))
151
152 exu.Print("[pOpt, vOpt]=", [pOpt, vOpt])
153 u = vOpt
154 exu.Print("optimum=",u)
155 exudynTestGlobals.testError = u - 0.0030262381366063158 #until 2022-02-20(changed to storeInternal): (0.0030262381385228617) #2020-12-18: (nElements=32) -2.7613614363986017e-05
156 exudynTestGlobals.testResult = u
157
158 if useGraphics and False:
159 # from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
160 import matplotlib.pyplot as plt
161
162 plt.close('all')
163 [figList, axList] = PlotOptimizationResults2D(pList, values, yLogScale=True)
164
165
166 #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
167 #ParameterVariation
168 [pList,v]=ParameterVariation(parameterFunction = ParameterFunction,
169 parameters = {'mass':(1,10,2), 'spring':(100,10000,3)},
170 useLogSpace=True,addComputationIndex=True,
171 showProgress=False)
172 #exu.Print("vList=", v)
173 u=v[3]
174 exudynTestGlobals.testError += u - 0.09814894553165972 #until 2022-02-20(changed to storeInternal):(0.09814894553377107) #2020-12-18: (nElements=32) -2.7613614363986017e-05
175 exudynTestGlobals.testResult += u
176 exu.Print('geneticOptimizationTest testResult=', exudynTestGlobals.testResult)
177 exu.Print('geneticOptimizationTest error=', exudynTestGlobals.testError)