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