ComputeSensitivitiesExample.py
You can view and download this file on Github: ComputeSensitivitiesExample.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: This example calculates the sensitivities of a mass-spring-damper-system
5# by varying mass, spring, ... and computing the forward/central difference
6# The output parameters used are the average absolute value of the displacement
7# and the static displacement
8#
9# Author: Peter Manzl, based on code from Johannes Gerstmayr
10# Date: 2022-02-10
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 ComputeSensitivities, PlotSensitivityResults
19from exudyn.utilities import AddSensorRecorder
20
21import numpy as np #for postprocessing
22
23#this is the function which is repeatedly called from ComputeSensitivities
24#parameterSet contains dictinary with varied parameters
25def ParameterFunction(parameterSet):
26 SC = exu.SystemContainer()
27 mbs = SC.AddSystem()
28
29
30 class P: pass #create 'namespace'
31
32 #default values
33 P.L=0.5 #spring length (for drawing)
34 P.mass = 1.6 #mass in kg
35 P.spring = 4000 #stiffness of spring-damper in N/m
36 P.damper = 8 #old: 8; damping constant in N/(m/s)
37 P.u0=-0.08 #initial displacement
38 P.v0=1 #initial velocity
39 P.force =80 #force applied to mass
40 P.computationIndex = 'Ref'
41
42 #update parameters:
43 for key in parameterSet: #includes empty dict!
44 setattr(P, key, parameterSet[key])
45
46
47 x0= P.force/P.spring #static displacement
48
49 #node for 3D mass point:
50 n1=mbs.AddNode(Point(referenceCoordinates = [P.L,0,0],
51 initialCoordinates = [P.u0,0,0],
52 initialVelocities= [P.v0,0,0]))
53
54 #ground node
55 nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
56
57 #add mass point (this is a 3D object with 3 coordinates):
58 massPoint = mbs.AddObject(MassPoint(physicsMass = P.mass, nodeNumber = n1))
59
60 #marker for ground (=fixed):
61 groundMarker=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
62 #marker for springDamper for first (x-)coordinate:
63 nodeMarker =mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
64
65 #spring-damper between two marker coordinates
66 nC = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundMarker, nodeMarker],
67 stiffness = P.spring, damping = P.damper))
68
69 #add load:
70 mbs.AddLoad(LoadCoordinate(markerNumber = nodeMarker,
71 load = P.force))
72 #add sensor:
73 fileName = 'solution/paramVarDisplacement'+ str(P.computationIndex) +'.txt'#
74 flagWriteFile = False
75 if P.computationIndex == 'Ref': flagWriteFile = True
76 sData = mbs.AddSensor(SensorObject(objectNumber=nC, fileName=fileName,
77 outputVariableType=exu.OutputVariableType.Force,
78 writeToFile=flagWriteFile))
79
80
81 steps = 1000 #number of steps to show solution
82 tEnd = 1 #end time of simulation
83
84 simulationSettings = exu.SimulationSettings()
85 simulationSettings.solutionSettings.writeSolutionToFile = False
86 simulationSettings.solutionSettings.sensorsWritePeriod = 5e-3 #output interval of sensors
87 simulationSettings.timeIntegration.numberOfSteps = steps
88 simulationSettings.timeIntegration.endTime = tEnd
89
90 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #no damping
91 if not(flagWriteFile):
92 sRecorder = AddSensorRecorder(mbs, sData, tEnd, simulationSettings.solutionSettings.sensorsWritePeriod, sensorOutputSize=1)
93
94 #exu.StartRenderer() #start graphics visualization
95 #mbs.WaitForUserToContinue() #wait for pressing SPACE bar to continue
96 mbs.Assemble()
97
98 #start solver:
99 mbs.SolveDynamic(simulationSettings)
100
101
102
103 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
104 #evaluate difference between reference and optimized solution
105 #reference solution:
106
107
108 if flagWriteFile:
109 data = np.loadtxt(fileName, comments='#', delimiter=',')
110 else:
111 data = mbs.variables['sensorRecord0']
112
113 avgPos = np.average(np.abs(data))
114 #+++++++++++++++++++++++++++++++++++++++++++++++++++++
115 #compute exact solution:
116 if False:
117 from matplotlib import plt
118
119 plt.close('all')
120 plt.plot(data[:,0], data[:,1], 'b-', label='displacement (m)')
121
122 ax=plt.gca() # get current axes
123 ax.grid(True, 'major', 'both')
124 ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
125 ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
126 plt.legend() #show labels as legend
127 plt.tight_layout()
128 plt.show()
129
130 return avgPos, x0
131
132
133
134#now perform the sensitivity analysis
135if __name__ == '__main__': #include this to enable parallel processing
136 import time
137 useMultiProcessing = exudyn
138 start_time = time.time()
139 n = [2, 2]
140 fVar = [1e-3, 1.5e-3, 1]
141 mRef = 1.5
142 kRef = 4000
143 [pList, valRef, valuesSorted, sensitivity] = ComputeSensitivities(parameterFunction=ParameterFunction,
144 parameters = {'mass': (mRef, fVar[0], n[0]),
145 'spring': (kRef,fVar[1], n[1]),
146 },
147 scaledByReference=False,
148 debugMode=True,
149 addComputationIndex=True,
150 useMultiProcessing=False,
151 showProgress=True,)
152
153 testResult = np.average(np.abs(sensitivity))
154 if True:
155 print("--- %s seconds ---" % (time.time() - start_time))
156 PlotSensitivityResults(valRef, valuesSorted, sensitivity, strYAxis=['avg. $|x|$', 'x0', ''])
157 else:
158 exu.Print('result of ConvexContactTest=',testResult)