minimizeExample.py
You can view and download this file on Github: minimizeExample.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: This example performs an 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# NOTE: using scipy.minimize with interface from Stefan Holzinger
9#
10# Author: Johannes Gerstmayr
11# Date: 2020-11-18
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.processing import Minimize, PlotOptimizationResults2D
20
21import numpy as np #for postprocessing
22import os
23from time import sleep
24
25
26#this is the function which is repeatedly called from ParameterVariation
27#parameterSet contains dictinary with varied parameters
28def ParameterFunction(parameterSet):
29 SC = exu.SystemContainer()
30 mbs = SC.AddSystem()
31
32 #default values
33 mass = 1.6 #mass in kg
34 spring = 4000 #stiffness of spring-damper in N/m
35 damper = 8 #old: 8; damping constant in N/(m/s)
36 u0=-0.08 #initial displacement
37 v0=1 #initial velocity
38 force =80 #force applied to mass
39
40 #process parameters
41 if 'mass' in parameterSet:
42 mass = parameterSet['mass']
43
44 if 'spring' in parameterSet:
45 spring = parameterSet['spring']
46
47 if 'force' in parameterSet:
48 force = parameterSet['force']
49
50 iCalc = 'Ref' #needed for parallel computation ==> output files are different for every computation
51 if 'computationIndex' in parameterSet:
52 iCalc = str(parameterSet['computationIndex'])
53 # print('iCAlc=', iCalc)
54
55
56 #mass-spring-damper system
57 L=0.5 #spring length (for drawing)
58
59 #node for 3D mass point:
60 n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0],
61 initialCoordinates = [u0,0,0],
62 initialVelocities= [v0,0,0]))
63
64 #ground node
65 nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
66
67 #add mass point (this is a 3D object with 3 coordinates):
68 massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n1))
69
70 #marker for ground (=fixed):
71 groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
72 #marker for springDamper for first (x-)coordinate:
73 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
74
75 #spring-damper between two marker coordinates
76 nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
77 stiffness = spring, damping = damper))
78
79 #add load:
80 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
81 load = force))
82 #add sensor:
83 sensorFileName = 'solution/paramVarDisplacement'+iCalc+'.txt'
84 mbs.AddSensor(SensorObject(objectNumber=nC, fileName=sensorFileName,
85 outputVariableType=exu.OutputVariableType.Displacement))
86 # print("sensorFileName",sensorFileName)
87
88 #print(mbs)
89 mbs.Assemble()
90
91 steps = 1000 #number of steps to show solution
92 tEnd = 1 #end time of simulation
93
94 simulationSettings = exu.SimulationSettings()
95 #simulationSettings.solutionSettings.solutionWritePeriod = 5e-3 #output interval general
96 simulationSettings.solutionSettings.writeSolutionToFile = False
97 simulationSettings.solutionSettings.sensorsWritePeriod = 2e-3 #output interval of sensors
98 simulationSettings.timeIntegration.numberOfSteps = steps
99 simulationSettings.timeIntegration.endTime = tEnd
100
101 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping
102
103 mbs.SolveDynamic(simulationSettings)
104
105 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
106 #evaluate difference between reference and optimized solution
107 #reference solution:
108 dataRef = np.loadtxt('solution/paramVarDisplacementRef.txt', comments='#', delimiter=',')
109 data = np.loadtxt(sensorFileName, comments='#', delimiter=',')
110
111 diff = data[:,1]-dataRef[:,1]
112
113 errorNorm = np.sqrt(np.dot(diff,diff))/steps*tEnd
114 #errorNorm = np.sum(abs(diff))/steps*tEnd
115
116 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
117 #draw solution (not during optimization!):
118 if 'plot' in parameterSet:
119
120 print('parameters=',parameterSet)
121 print('file=', sensorFileName)
122 print('error=', errorNorm)
123 import matplotlib.pyplot as plt
124 from matplotlib import ticker
125
126 plt.close('all')
127 plt.plot(dataRef[:,0], dataRef[:,1], 'b-', label='Ref, u (m)')
128 plt.plot(data[:,0], data[:,1], 'r-', label='u (m)')
129
130 ax=plt.gca() # get current axes
131 ax.grid(True, 'major', 'both')
132 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
133 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
134 plt.legend() #show labels as legend
135 plt.tight_layout()
136 plt.show()
137
138
139 if True: #not needed in serial version
140 if iCalc != 'Ref':
141 os.remove(sensorFileName) #remove files in order to clean up
142 while(os.path.exists(sensorFileName)): #wait until file is really deleted -> usually some delay
143 sleep(0.001) #not nice, but there is no other way than that
144
145 del mbs
146 del SC
147
148 # print(parameterSet, errorNorm)
149 return errorNorm
150
151
152#now perform parameter variation
153if __name__ == '__main__': #include this to enable parallel processing
154 import time
155
156 refval = ParameterFunction({}) # compute reference solution
157 print("refval =", refval)
158 if False:
159 #val2 = ParameterFunction({'mass':1.6, 'spring':4000, 'force':80, 'computationIndex':0, 'plot':''}) # compute reference solution
160 val2 = ParameterFunction({'mass': 1.7022816582583309, 'spring': 4244.882757974497, 'force': 82.62761337061548, 'computationIndex':0, 'plot':''}) # compute reference solution
161 #val2 = ParameterFunction({, 'computationIndex':0, 'plot':''}) # compute reference solution
162
163 if True:
164 #the following settings give approx. 6 digits accuraet results after 167 iterations
165 start_time = time.time()
166 [pOpt, vOpt, pList, values] = Minimize(objectiveFunction = ParameterFunction,
167 parameters = {'mass':(1,10), 'spring':(100,10000), 'force':(1,250)}, #parameters provide search range
168 showProgress=True,
169 debugMode=False,
170 addComputationIndex = True,
171 tol = 1e-1, #this is a internal parameter, not directly coupled loss
172 options={'maxiter':200},
173 resultsFile='solution/test.txt'
174 )
175 print("--- %s seconds ---" % (time.time() - start_time))
176
177 print("optimum parameters=", pOpt)
178 print("minimum value=", vOpt)
179
180 from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import
181 import matplotlib.pyplot as plt
182 #from matplotlib import cm
183 #from matplotlib.ticker import LinearLocator, FormatStrFormatter
184 import numpy as np
185 colorMap = plt.cm.get_cmap('jet') #finite element colors
186
187 #for negative values:
188 if min(values) <= 0:
189 values = np.array(values)-min(values)*1.001+1e-10
190
191 plt.close('all')
192 [figList, axList] = PlotOptimizationResults2D(pList, values, yLogScale=True)