black-box RC model

PhysicalModels/genRCData.m

@@ -11,7 +11,7 @@ duration = 10;
 resolution = 0.001;
 N = duration / resolution;
 t = 0:resolution:(duration-resolution);
-noise = 250; % signal-to-noise ratio in dB
+noise = 25; % signal-to-noise ratio in dB
 peak = 10;  % amplitude of signal
 freq = 1;   % frequency of sine signal
 
@@ -20,19 +20,22 @@ in = t' * peak / duration;
 % sim_in = [t', awgn(in, noise, 'measured')];
 sim_in = [t', in];
 data = runSim('RC', sim_in);
-save('resRC_ramp.mat', 'data', 'peak', 'resolution');
+save("resRC_ramp_" + num2str(noise) + ".mat",...
+    'data', 'peak', 'resolution');
 
 % sine input
 in = sin(2*pi*freq*t') * peak;
 % sim_in = [t', awgn(in, noise, 'measured')];
 sim_in = [t', in];
 data = runSim('RC', sim_in);
-save('resRC_sine.mat', 'data', 'freq', 'peak', 'resolution');
+save("resRC_sin_" + num2str(noise) + ".mat",...
+    'data', 'freq', 'peak', 'resolution');
 
 % combined ramp + sine + noise
 in = (t' * peak / duration) + sin(2*pi*freq*t') * peak;
 sim_in = [t', awgn(in, noise, 'measured')];
 data = runSim('RC', sim_in);
-save('resRC_comb.mat', 'data', 'freq', 'peak', 'resolution');
+save("resRC_comb_" + num2str(noise) +".mat",...
+    'data', 'freq', 'peak', 'resolution');
 
 clear data freq peak N t duration resolution in sim_in
\ No newline at end of file

blackBoxRC.m

+% Estimates a LSTM model of a RC circuit.
+% 
+% See:
+% - https://www.mathworks.com/help/deeplearning/ug/long-short-term-memory-networks.html
+% - https://www.mathworks.com/help/deeplearning/examples/sequence-to-sequence-regression-using-deep-learning.html
+% - https://www.mathworks.com/help/deeplearning/examples/time-series-forecasting-using-deep-learning.html
+
+close all;
+
+% Load training data
+rcdat = load('resRC_comb_250.mat');
+t = rcdat.data(:,1);    % timestamps
+u = rcdat.data(:,2);    % input
+y = rcdat.data(:,3) * 1000;    % output
+
+mu = mean(u);
+sig = std(u);
+stdU = (u - mu) / sig;
+
+% Split into sub-sequences
+seqLength = 500;
+seqU = reshape(stdU, [], seqLength);
+seqY = reshape(y, [], seqLength);
+
+
+XTrain = mat2cell(seqU, ones(1, size(seqU,1)));
+YTrain = mat2cell(seqY, ones(1, size(seqY,1)));
+
+% creating LSTM architecture
+numFeatures = 1;
+numHiddenUnits = 16;
+numResponses = 1;
+
+layers = [ ...
+    sequenceInputLayer(numFeatures)
+    lstmLayer(numHiddenUnits,'OutputMode','sequence')
+    lstmLayer(numHiddenUnits,'OutputMode','sequence')
+    lstmLayer(numHiddenUnits,'OutputMode','sequence')
+    lstmLayer(numHiddenUnits,'OutputMode','sequence')
+    fullyConnectedLayer(numResponses)
+    regressionLayer];
+
+options = trainingOptions('adam', ...
+    'MaxEpochs',40, ...
+    'MiniBatchSize', 5, ...
+    'GradientThreshold',1, ...
+    'InitialLearnRate',0.1, ...
+    'LearnRateSchedule','piecewise', ...
+    'LearnRateDropPeriod',5, ...
+    'LearnRateDropFactor',0.7, ...
+    'Verbose',0, ...
+    'Plots','training-progress');
+
+net = trainNetwork(XTrain, YTrain, layers, options);
+
+YPred = cell2mat(predict(net, {u'}));
+plot(t, y, t, YPred);
+legend('Actual', 'Prediction');
\ No newline at end of file

docs/0_background.md

@@ -15,10 +15,15 @@ Grey-box models add some obfuscation to the system dynamics. While the structure
 
 An example could be a simple RC circuit with unknown capacitance and resistance. The equations governing the current through the components are known. A physical model can be constructed based off of those equations. The model, in combination with actual measurements from the circuit then estimates the values of parameters in the equations such that the difference between the model's predictions and the actual readings is minimized.
 
+### State-space equations as grey-boxes
+
 ## Black-box models
 
 A Black-box model foregoes any *apriori* knowledge about the distribution of system parameters. Instead it learns the mechanics from scratch. A neural network used to approximate an RC circuit is a black-box model. The network simply learns the mappings from the inputs to the outputs.
 
+### LSTM networks as black-boxes
+
+Long short-term memory networks are a variant of recurrent neural networks. They are *stateful* models. An LSTM network is able to memorize earlier states such that it affects later states. This is important when modelling dynamic systems wherre the future states of the system depend on the rate of change of states or other higher order phenomenon.
 
 ### Further reading:
 

docs/1_RL_Circuit.md

@@ -77,5 +77,9 @@ We set `R` as the learnable parameter under various conditions (noisy inputs, wi
 
 The results show that the model is highly susceptible to noise. The mathematical structure of the model does not account for a parasitic element, therefore the peak accuracy it is able to achieve theoretically tops out before reaching a 100%. However, the limited model still arrives at the correct estimate for `R`. The worst case for this grey-box model is under noisy training data and unaccounted-for electrical effects which reduce the model performance to 0.
 
+### Black box model
+
+
+
 [1]: https://en.wikibooks.org/wiki/Control_Systems#Modern_Control_Methods
 [2]: https://en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws#Kirchhoff's_voltage_law_(KVL)
\ No newline at end of file

dummyLSTMone.m

+% Training an LSTM network to predict the next number
+% from an input sequence.
+% 
+% Input:  1,2,3
+% Output: 4 
+N = 1000;
+NTest = 200;
+seqLen = 10;
+
+seq = reshape([1:N], seqLen, [])';
+testSeq = reshape([N+1:N+NTest], seqLen, [])';
+
+% The input variables are normalized to have 0 mean
+% and 1 variance.
+mu = mean(seq(:));
+sig = std(seq(:));
+stdSeq = (seq - mu) / sig;
+stdTest = (testSeq - mu) / sig;
+
+% Preparing input data as a cell array of row vectors.
+XTrain = mat2cell(stdSeq(1:end-1,:),...
+                  ones(1, N/seqLen - 1));
+XTest = mat2cell(stdTest(1:end-1,:),...
+                  ones(1, NTest/seqLen - 1));
+
+% For format of the response, see
+% https://www.mathworks.com/help/deeplearning/ref/trainnetwork.html
+% The labels are a matrix of the next number for each
+% training sequence.
+YTrain = seq(2:end,1);
+YTest = testSeq(2:end,1);
+
+% Defining an LSTM network. An LSTM layer takes as
+% input `numHiddenUnits` which is the number of time
+% steps that an LSTM cell "memorizes" i.e. the window
+% over which the LSTM layer operates.
+layers = [...
+    sequenceInputLayer(1)
+    lstmLayer(10,'OutputMode','last')
+    fullyConnectedLayer(1)
+    regressionLayer];
+
+options = trainingOptions('adam', ...
+    'MaxEpochs',150, ...
+    'MiniBatchSize', 5, ...
+    'GradientThreshold',1, ...
+    'InitialLearnRate',1.0, ...
+    'LearnRateSchedule','piecewise', ...
+    'LearnRateDropPeriod', 20, ...
+    'LearnRateDropFactor',0.5, ...
+    'Verbose',0, ...
+    'Plots','training-progress');
+
+net = trainNetwork(XTrain, YTrain, layers, options);
+
+YPred = int32(predict(net, XTrain));
\ No newline at end of file

dummyLSTMseq.m

+% Training an LSTM network to predict a sequence
+% of numbers following an input sequence.
+% 
+% Input:  1,2,3
+% Output: 2,3,4
+N = 1000;
+NTest = 200;
+seqLen = 10;
+
+seq = reshape([1:N], seqLen, [])';
+testSeq = reshape([N+1:N+NTest], seqLen, [])';
+
+% The input variables are normalized to have 0 mean
+% and 1 variance.
+mu = mean(seq(:));
+sig = std(seq(:));
+stdSeq = (seq - mu) / sig;
+stdTest = (testSeq - mu) / sig;
+
+% Preparing input data as a cell array of row vectors.
+XTrain = mat2cell(stdSeq,ones(1, N/seqLen));
+XTest = mat2cell(stdTest, ones(1, NTest/seqLen));
+
+% For format of the response, see
+% https://www.mathworks.com/help/deeplearning/ref/trainnetwork.html
+% The labels are a cell array of sequences offset by 1
+YTrain = mat2cell(seq + 1, ones(1, N/seqLen));
+YTest = mat2cell(testSeq + 1, ones(1, NTest/seqLen));
+
+% Defining an LSTM network. An LSTM layer takes as
+% input `numHiddenUnits` which is the number of time
+% steps that an LSTM cell "memorizes" i.e. the window
+% over which the LSTM layer operates.
+layers = [...
+    sequenceInputLayer(1)
+    lstmLayer(10,'OutputMode','sequence')
+    fullyConnectedLayer(1)
+    regressionLayer];
+
+options = trainingOptions('adam', ...
+    'MaxEpochs',150, ...
+    'MiniBatchSize', 5, ...
+    'GradientThreshold',1, ...
+    'InitialLearnRate',1.0, ...
+    'LearnRateSchedule','piecewise', ...
+    'LearnRateDropPeriod',20, ...
+    'LearnRateDropFactor',0.5, ...
+    'Verbose',0, ...
+    'Plots','training-progress');
+
+net = trainNetwork(XTrain, YTrain, layers, options);
+
+YPred = int32(cell2mat(predict(net, XTrain)));
\ No newline at end of file

greyBoxRC.m

@@ -21,7 +21,7 @@ data = iddata(y, u, rcdat.resolution);
 % Run a grey box estimate of the system using data
 sys = greyest(data, sys_init);
 
-fprintf('Parasysmeter value: %.5f\n', getpvec(sys));
+fprintf('Parameter value: %.5f\n', getpvec(sys));
 
 % state-space representation of RC circuit
 function [A, B, C, D] = ssRC(R, c, ~)