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