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