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