tutorialNeuralNetwork.py
You can view and download this file on Github: tutorialNeuralNetwork.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: tutorial for machine learning with Exudyn;
5# correct positioning of rigid body mounted on strings by model trained with machine learning
6# data is created with static computations, then inverse model is trained with pytorch
7#
8# Model: A rigid body with height 0.4, width 0.2 and depth 0.2 and 1 kg, mounted on two soft strings with stiffness 500N/m
9#
10# Author: Johannes Gerstmayr
11# Date: 2024-02-16
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.utilities import *
19# from exudyn.signalProcessing import GetInterpolatedSignalValue
20
21import sys
22import numpy as np
23# #from math import sin, cos, sqrt,pi
24import os
25os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE" #for multiprocessing problems with pytorch
26
27from numpy.random import rand
28np.random.seed(0)
29
30doTraining = True
31
32#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++
33#create an Exudyn model of a rigid body on two strings
34
35SC = exu.SystemContainer()
36mbs = SC.AddSystem()
37
38
39np.random.seed(0)
40
41#Geometry:
42#[0,0,0] is at bottom of left tower
43#[L,0,0] is at bottom of right tower
44#[L,H,0] is at top of left tower
45
46L = 4 #m; distance of columns
47H = 3 #m; height of columns
48w = 0.2 #m; size of rigid body
49m = 1 #kg; weight of mass
50g = 9.81 #m/s^2; gravity
51
52sideLengths = [w,w*2,w]
53density = m/w**3 #average density of rigid body with dimensions
54
55#left and right tower
56pTower0 = np.array([0,H,0])
57pTower1 = np.array([L,H,0])
58#local attachment points at mass
59localPosMass0=[-0.5*w,w,0]
60localPosMass1=[ 0.5*w,w,0]
61#center location of rigid body
62#pRigidMid = np.array([0.5*L, 0.5*H,0])
63pRigidMid = np.array([0.5*L, 0.5*H,0])
64# pInit = [1.2,0.78,0]
65pInit = pRigidMid
66
67k = 500 #string stiffness
68d = 0.01*k #string damping for dynamic simulation
69tEnd = 1
70stepSize = 0.01
71
72gGround = [GraphicsDataOrthoCubePoint(centerPoint=[0,0.5*H,0],size=[0.5*w,H,w],color=color4grey)]
73gGround += [GraphicsDataOrthoCubePoint(centerPoint=[L,0.5*H,0],size=[0.5*w,H,w],color=color4grey)]
74gGround += [GraphicsDataOrthoCubePoint(centerPoint=[0.5*L,-0.5*w,0],size=[L,w,w],color=color4grey)]
75oGround = mbs.CreateGround(graphicsDataList=gGround)
76
77gRigid = [GraphicsDataOrthoCubePoint(size=sideLengths, color=color4dodgerblue)]
78oRigid = mbs.CreateRigidBody(referencePosition=pInit,
79 inertia = InertiaCuboid(density, sideLengths),
80 gravity = [0,-g,0],
81 graphicsDataList=gRigid)
82nRigid = mbs.GetObject(oRigid)['nodeNumber'] #used later
83
84oString0 = mbs.CreateSpringDamper(bodyList=[oGround, oRigid],
85 localPosition0=[0,H,0],
86 localPosition1=localPosMass0,
87 stiffness = k,
88 damping = d,
89 drawSize = 0, #draw as line
90 )
91oString1 = mbs.CreateSpringDamper(bodyList=[oGround, oRigid],
92 localPosition0=[L,H,0],
93 localPosition1=localPosMass1,
94 stiffness = k,
95 damping = d,
96 drawSize = 0, #draw as line
97 )
98sRigid = mbs.AddSensor(SensorBody(bodyNumber=oRigid, storeInternal=True,
99 outputVariableType=exu.OutputVariableType.Position))
100
101# compute string lengths for given rigid body center position in straight configuration
102# used for initialization of static computation
103def ComputeStringLengths(pRigid):
104 L0 = np.array(pRigid)+localPosMass0-pTower0
105 L1 = np.array(pRigid)+localPosMass1-pTower1
106
107 return [NormL2(L0), NormL2(L1)]
108
109
110mbs.Assemble()
111simulationSettings = exu.SimulationSettings() #takes currently set values or default values
112simulationSettings.solutionSettings.writeSolutionToFile = not doTraining
113
114# # this leads to flipped results => good example !
115# simulationSettings.staticSolver.numberOfLoadSteps = 10
116# simulationSettings.staticSolver.stabilizerODE2term = 2 #add virtual stiffness due to mass; helps static solver to converge for such cases
117
118simulationSettings.staticSolver.numberOfLoadSteps = 2
119simulationSettings.staticSolver.stabilizerODE2term = 20 #add virtual stiffness due to mass; helps static solver to converge for such cases
120simulationSettings.staticSolver.computeLoadsJacobian = False #due to bug in loadsJacobian
121simulationSettings.staticSolver.verboseMode = 0
122
123if False:
124 tEnd = 20 #for visualization of dynamic case
125 simulationSettings.solutionSettings.sensorsWritePeriod = stepSize
126 # simulationSettings.timeIntegration.simulateInRealtime = True
127 simulationSettings.timeIntegration.endTime = tEnd
128 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
129
130 exu.StartRenderer()
131 mbs.WaitForUserToContinue()
132
133 mbs.SolveStatic(simulationSettings)
134 SC.WaitForRenderEngineStopFlag()
135 mbs.SolveDynamic(simulationSettings)
136
137 SC.WaitForRenderEngineStopFlag()
138 exu.StopRenderer()
139
140 mbs.PlotSensor(sRigid, components=[0,1,2]) #plot vertical displacement
141
142 #this shows the deviation due to string stiffness and rotation of rigid body
143 print('final pos=',mbs.GetSensorValues(sRigid))
144
145 sys.exit()
146
147#this function is called to compute real position from ideal position p
148def ComputePositionFromStringLengths(p):
149 [L0,L1] = ComputeStringLengths(p)
150 refCoordsRigid[0:3] = p #override position
151 mbs.SetNodeParameter(nodeNumber=nRigid,
152 parameterName='referenceCoordinates',
153 value=refCoordsRigid)
154 mbs.SetObjectParameter(objectNumber=oString0,
155 parameterName='referenceLength',
156 value=L0)
157 mbs.SetObjectParameter(objectNumber=oString1,
158 parameterName='referenceLength',
159 value=L1)
160 mbs.Assemble()
161
162 try:
163 mbs.SolveStatic(simulationSettings)
164 positionList.append(p) #this is the ideal position; used to calculate deviation
165 #we map targeted (real) positions to original lengths
166 realPos = mbs.GetSensorValues(sRigid)
167 #compute original lengths to realPos
168 [L0orig,L1orig] = ComputeStringLengths(realPos)
169
170 #targetList.append([L0,L1]) #we only need to store the deviation
171 diff = [L0-L0orig,L1-L1orig]
172 return [realPos, diff, [L0, L1]]
173 except:
174 print('solver failed for:',p,',Ls=',[L0,L1])
175 return [None]
176
177#%%+++++++++++++++++++++++++++++++++++++++++++++
178#now create data for desired positions
179gridX = 20*2
180gridY = 20*2
181nExamples = gridX*gridY
182nExamplesTest = int(nExamples*0.1)
183pRangeX = 2.4
184pRangeY = 2.4
185positionList = []
186inputList = []
187targetList = []
188#store reference coordinates for rotations
189refCoordsRigid = mbs.GetNodeParameter(nodeNumber=nRigid, parameterName='referenceCoordinates')
190
191gridValues=np.zeros((gridX,gridY,4))
192i=0
193ix = 0
194iy = 0
195while i < nExamples+nExamplesTest:
196 if i < nExamples:
197 ix = i%gridX
198 iy = int(i/gridX)
199 x0 = pRangeX*(ix/gridX-0.5)
200 y0 = pRangeY*(iy/gridY-0.5)
201 else:
202 x0 = pRangeX*(rand()-0.5)
203 y0 = pRangeY*(rand()-0.5)
204
205 p = pRigidMid + [x0, y0, 0]
206
207 rv = ComputePositionFromStringLengths(p)
208
209 if rv[0] is not None:
210 realPos = rv[0]
211 diff = rv[1]
212 [L0,L1] = rv[2]
213 targetList.append(diff)
214 inputList.append(list(realPos)[0:2]) #correction on position
215
216 if i < nExamples:
217 gridValues[ix,iy,0:2] = diff
218 gridValues[ix,iy,2:4] = realPos[0:2]
219
220 if max(diff)>0.5:
221 print('++++++++++++++++++++++++++')
222 print('ideal pos=',p)
223 print('realPos =',realPos)
224 print('lengths =',L0,L1)
225
226 i += 1
227
228 #else: do not increment and retry ...
229
230
231inputsExudynList = inputList[0:nExamples]
232targetsExudynList = targetList[0:nExamples]
233inputsExudynTestList = inputList[nExamples:]
234targetsExudynTestList = targetList[nExamples:]
235
236print('created',nExamples+nExamplesTest,'samples')
237
238#%%
239def PlotGridFunction(data, name='error', closeAll=False):
240 import numpy as np
241 import matplotlib.pyplot as plt
242 from mpl_toolkits.mplot3d import Axes3D
243
244 if closeAll: plt.close('all')
245
246 # Generate grid coordinates (if not already generated)
247 # For example, if gridX and gridY are the dimensions of your grid
248 gridX, gridY = data.shape
249 x = np.linspace( -0.5*pRangeX, 0.5*pRangeX, gridX)
250 y = np.linspace(-0.5*pRangeY, 0.5*pRangeY, gridY)
251 X, Y = np.meshgrid(x, y)
252
253 # Create the contour plot
254 fig=plt.figure(figsize=(8, 6))
255 ax = fig.add_subplot(111, projection='3d')
256
257 # contour = plt.contourf(X, Y, dataX, cmap='hot', levels=100)
258 contour = ax.plot_surface(X, Y, data, cmap='viridis')
259
260 plt.colorbar(contour)
261 plt.title(name)
262 plt.xlabel('X-axis')
263 plt.ylabel('Y-axis')
264
265 # Display the plot
266 plt.show()
267
268PlotGridFunction(gridValues[:, :, 0], name='error X', closeAll=True)
269PlotGridFunction(gridValues[:, :, 1], name='error Y')
270
271
272if not doTraining:
273 sys.exit()
274
275
276
277
278
279#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
280#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
281
282#now train
283import torch
284import torch.nn as nn
285import torch.optim as optim
286from torch.utils.data import DataLoader, TensorDataset
287
288hiddenLayerSize = 8*4 # Example size, adjust as needed
289batchSize = 16*8
290learningRate = 0.002
291nTrainingEpochs = 1000*10
292lossThreshold = 0.0002
293
294# torch.set_num_threads(1)
295#prepare data:
296
297# Convert lists to PyTorch tensors
298dtype=torch.float32
299inputsExudyn = torch.tensor(inputsExudynList, dtype=dtype)
300targetsExudyn = torch.tensor(targetsExudynList, dtype=dtype)
301inputsExudynTest = torch.tensor(inputsExudynTestList, dtype=dtype)
302targetsExudynTest = torch.tensor(targetsExudynTestList, dtype=dtype)
303
304# Create TensorDatasets
305train_dataset = TensorDataset(inputsExudyn, targetsExudyn)
306test_dataset = TensorDataset(inputsExudynTest, targetsExudynTest)
307
308# Create DataLoaders
309trainLoader = DataLoader(train_dataset, batch_size=batchSize, shuffle=True)
310testLoader = DataLoader(test_dataset, batch_size=batchSize, shuffle=False)
311
312class MultiplicativeLayer(nn.Module):
313 def __init__(self):
314 super(MultiplicativeLayer, self).__init__()
315 # You can also include parameters or learnable weights here if necessary
316
317 def forward(self, x):
318 # Assuming x is of shape [batch_size, num_features]
319 # We'll create new features by multiplying pairs of features
320 new_features = []
321 for i in range(x.size(1)):
322 for j in range(i + 1, x.size(1)):
323 new_features.append(x[:, i] * x[:, j])
324
325 # Concatenate the original features with the new, multiplied features
326 if new_features:
327 new_features = torch.stack(new_features, dim=1)
328 x = torch.cat((x, new_features), dim=1)
329 return x
330
331class ModelNN(nn.Module):
332 def __init__(self, input_size, hiddenLayerSize, output_size):
333 super(ModelNN, self).__init__()
334 self.multiplicativeLayer = MultiplicativeLayer()
335 # self.fc1 = nn.Linear(input_size + (input_size * (input_size - 1) // 2),
336 # hiddenLayerSize,dtype=dtype) # Adjusted input size
337 self.fc1 = nn.Linear(input_size, hiddenLayerSize,dtype=dtype)
338 self.relu = nn.ReLU()
339 self.leakyRelu= nn.LeakyReLU()
340 self.elu = nn.ELU()
341 # self.relu = nn.Sigmoid()
342 # self.relu = nn.Tanh()
343 self.fc2 = nn.Linear(hiddenLayerSize, hiddenLayerSize,dtype=dtype)
344 self.fc3 = nn.Linear(hiddenLayerSize, hiddenLayerSize,dtype=dtype)
345 self.lastLayer = nn.Linear(hiddenLayerSize, output_size,dtype=dtype)
346
347 def forward(self, x):
348 # x = self.multiplicativeLayer(x) #not a great help
349 x = self.fc1(x)
350 x = self.relu(x)
351 x = self.fc2(x)
352 x = self.leakyRelu(x)
353 x = self.fc3(x)
354 x = self.elu(x)
355 x = self.lastLayer(x)
356 return x
357
358
359
360input_size = inputsExudyn.shape[1] # Number of input features
361output_size = targetsExudyn.shape[1] # Assuming regression task, adjust for classification
362
363model = ModelNN(input_size, hiddenLayerSize, output_size)
364lossFunction = nn.MSELoss() # Mean Squared Error Loss for regression, adjust for classification
365optimizer = optim.Adam(model.parameters(), lr=learningRate,)
366
367#++++++++++++++++++++++++++++++++++++++++++++++++++++
368#set up training
369
370lossHistory = []
371minLoss = 1e10
372# Train the network
373for epoch in range(nTrainingEpochs): # 100 epochs
374 for inputs, targets in trainLoader:
375 optimizer.zero_grad()
376
377 # Forward pass
378 outputs = model(inputs)
379 loss = lossFunction(outputs, targets)
380
381 # Backward pass and optimization
382 loss.backward()
383 optimizer.step()
384
385 lossHistory.append([epoch, np.sqrt(loss.item())])
386 minLoss = min(minLoss, np.sqrt(loss.item()))
387 lossVal = np.sqrt(loss.item())
388
389 if lossVal < lossThreshold:
390 print(f'loss threshold reached at: epoch {epoch+1}/{nTrainingEpochs}, Loss: {lossVal}')
391 break
392
393 if (epoch+1) % 50 == 0:
394 print(f'Epoch {epoch+1}/{nTrainingEpochs}, Loss: {lossVal}')
395
396print('min loss=',minLoss)
397
398#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++
399#evaluate:
400model.eval() # Set the model to evaluation mode
401totalLoss = 0
402count = 0
403
404with torch.no_grad(): # No need to track gradients for evaluation
405 for inputs, targets in testLoader:
406 outputs = model(inputs)
407 loss = torch.sqrt(((outputs - targets) ** 2).mean()) # Calculating RMSE for each batch
408 totalLoss += loss.item()
409 count += 1
410
411averageRMSE = totalLoss / count
412
413# Call the evaluate_model function with the test_loader and your model
414print(f"\nTest RMSE: {averageRMSE}\n")
415
416for test in range(10):
417 x = inputsExudynTest[test:test+1]
418 print('x=',x,', shape',x.shape)
419 y = model(x).tolist()[0] #convert output to list
420
421 yRef = targetsExudynTest[test:test+1]
422
423 print('++++++++++++++++++++')
424 print('input: ',x,'\ntarget: ',yRef,'\npredict:',y)
425 # mbs.PlotSensor(result, components=[0], closeAll=test==0, newFigure=False,
426 # labels=['RNN'], yLabel='displacement (m)',
427 # colorCodeOffset=test)
428 # mbs.PlotSensor(result, components=[1], newFigure=False,
429 # labels=['reference'], yLabel='displacement (m)',
430 # colorCodeOffset=test,lineStyles=[':'])
431
432testErrorGrid = np.zeros((gridX, gridY))
433maxError = 0
434for iy in range(gridY):
435 for ix in range(gridX):
436 x0 = pRangeX*((ix+0.5)/gridX-0.5)*0.5
437 y0 = pRangeY*((iy+0.5)/gridY-0.5)*0.5
438
439 p = pRigidMid + [x0, y0, 0]
440
441 x = torch.tensor([list(p[0:2])], dtype=dtype)
442 y = model(x).tolist()[0]
443
444 rv = ComputePositionFromStringLengths(p)
445
446 diff = rv[1]
447 err = np.array(diff) - y
448 maxError = max(maxError, abs(err[0]))
449
450 testErrorGrid[ix,iy] = err[0]
451
452print('maxError', maxError)
453PlotGridFunction(testErrorGrid, name='test error X', closeAll=True)
454#PlotGridFunction(gridValues[:, :, 0], name='error X', closeAll=True)
455
456#%%+++++++++++++++
457if False:
458 lossData = np.array(lossHistory)
459 import matplotlib.pyplot as plt
460 from exudyn.plot import PlotSensor
461 PlotSensor(None,lossData,components=[0], closeAll=True, logScaleY=True)