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