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
 24testIsActive = exu.sys.get('testIsActive', False)
 25
 26dataRef = None
 27#this is the function which is repeatedly called from ParameterVariation
 28#parameterSet contains dictinary with varied parameters
 29def ParameterFunction(parameterSet):
 30    SC = exu.SystemContainer()
 31    mbs = SC.AddSystem()
 32    global dataRef #can be avoided, if reference solution is written to file!
 33
 34    
 35    #++++++++++++++++++++++++++++++++++++++++++++++
 36    #++++++++++++++++++++++++++++++++++++++++++++++
 37    #store default parameters in structure (all these parameters can be varied!)
 38    class P: pass #create emtpy structure for parameters; simplifies way to update parameters
 39
 40    P.mass = 1.6          #mass in kg
 41    P.spring = 4000       #stiffness of spring-damper in N/m
 42    P.damper = 8    #old: 8; damping constant in N/(m/s)
 43    P.u0=-0.08            #initial displacement
 44    P.v0=1                #initial velocity
 45    P.force =80           #force applied to mass
 46    P.computationIndex = ''
 47
 48    #now update parameters with parameterSet (will work with any parameters in structure P)
 49    for key,value in parameterSet.items():
 50        setattr(P,key,value)
 51
 52    #++++++++++++++++++++++++++++++++++++++++++++++
 53    #++++++++++++++++++++++++++++++++++++++++++++++
 54    #create model using parameters P starting here:
 55
 56    L=0.5               #spring length (for drawing)
 57    n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0], 
 58                         initialCoordinates = [P.u0,0,0], 
 59                         initialVelocities= [P.v0,0,0]))
 60    
 61    #ground node
 62    nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 63    
 64    #add mass point (this is a 3D object with 3 coordinates):
 65    massPoint = mbs.AddObject(MassPoint(mass = P.mass, nodeNumber = n1))
 66    
 67    #marker for ground (=fixed):
 68    groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 69    #marker for springDamper for first (x-)coordinate:
 70    nodeMarker  =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 71    
 72    #spring-damper between two marker coordinates
 73    nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker], 
 74                                              stiffness = P.spring, damping = P.damper)) 
 75    
 76    #add load:
 77    mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker, load = P.force))
 78    #add sensor:
 79    
 80    sDisp = mbs.AddSensor(SensorObject(objectNumber=nC, 
 81                               storeInternal=True,
 82                               outputVariableType=exu.OutputVariableType.Displacement))
 83 
 84    #++++++++++++++++++++++++++++++++++++++++++++++
 85    #assemble and compute solution
 86    mbs.Assemble()
 87    
 88    steps = 100  #number of steps to show solution
 89    tEnd = 1     #end time of simulation
 90    
 91    simulationSettings = exu.SimulationSettings()
 92    simulationSettings.solution.file.write = False
 93    simulationSettings.solution.sensors.writePeriod = 2e-3  #output interval of sensors
 94    simulationSettings.timeIntegration.numberOfSteps = steps
 95    simulationSettings.timeIntegration.endTime = tEnd
 96    
 97    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping
 98    
 99    mbs.SolveDynamic(simulationSettings)
100    
101    #+++++++++++++++++++++++++++++++++++++++++++++++++++++
102    #evaluate difference between reference and optimized solution
103    #reference solution:
104    data = mbs.GetSensorStoredData(sDisp)
105
106    if P.computationIndex=='Ref': 
107        dataRef = data
108        
109    diff = data[:,1]-dataRef[:,1]
110    
111    errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd
112
113    return errorNorm
114
115#generate reference data, also in multi-threaded version this will be computed for all threads!
116refval = ParameterFunction({'computationIndex':'Ref'}) # compute reference solution
117
118#now perform parameter variation
119if __name__ == '__main__': #include this to enable parallel processing
120    import time
121
122    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
123    #GeneticOptimization    
124    start_time = time.time()
125    [pOpt, vOpt, pList, values] = GeneticOptimization(objectiveFunction = ParameterFunction, 
126                                         parameters = {'mass':(1,10), 'spring':(100,10000), 'force':(1,1000)}, #parameters provide search range
127                                         numberOfGenerations = 2,
128                                         populationSize = 10,
129                                         elitistRatio = 0.1,
130                                         crossoverProbability = 0.1,
131                                         rangeReductionFactor = 0.7,
132                                         addComputationIndex=True,
133                                         randomizerInitialization=0, #for reproducible results
134                                         distanceFactor = 0.1, #for this example only one significant minimum
135                                         debugMode=False,
136                                         useMultiProcessing=False, #may be problematic for test
137                                         showProgress=False,
138                                         resultsFile = 'solution/geneticOptimizationTest.txt',
139                                         )
140    exu.Print("--- %s seconds ---" % (time.time() - start_time))
141
142    exu.Print("[pOpt, vOpt]=", [pOpt, vOpt])
143    u = vOpt
144    exu.Print("optimum=",u)
145    exu.sys['testResult'] = u
146
147    if not testIsActive and False:
148        # from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
149        import matplotlib.pyplot as plt
150
151        plt.close('all')
152        [figList, axList] = PlotOptimizationResults2D(pList, values, yLogScale=True)
153        
154
155    #%%++++++++++++++++++++++++++++++++++++++++++++++++++++
156    #ParameterVariation    
157    [pList,v]=ParameterVariation(parameterFunction = ParameterFunction, 
158                       parameters = {'mass':(1,10,2), 'spring':(100,10000,3)},
159                       useLogSpace=True,addComputationIndex=True,
160                       showProgress=False)
161    #exu.Print("vList=", v)
162    u=v[3]
163    exu.sys['testResult'] += u
164    exu.Print('geneticOptimizationTest testResult=', exu.sys['testResult'])