grey box estimation for RC circuit

PhysicalModels/genRCData.m

@@ -7,29 +7,32 @@
 % Other variables like peak, freq etc are stored for
 % relevant inputs.
 
-duration = 1000;
-resolution = 0.1;
+duration = 10;
+resolution = 0.001;
 N = duration / resolution;
 t = 0:resolution:(duration-resolution);
-noise = 25; % signal-to-noise ratio in dB
+noise = 250; % signal-to-noise ratio in dB
 peak = 10;  % amplitude of signal
 freq = 1;   % frequency of sine signal
-bias = 3;   % bias of offset sine signal
 
-% constant input
-in = ones(N,1) * peak;
-sim_in = [t', awgn(in, noise, 'measured')];
-datRC = runSim('RC', sim_in);
-save('resRC_const.mat', 'datRC', 'peak', 'noise');
+% ramp input
+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');
 
 % sine input
 in = sin(2*pi*freq*t') * peak;
-sim_in = [t', awgn(in, noise, 'measured')];
-datRC = runSim('RC', sim_in);
-save('resRC_sine.mat', 'datRC', 'freq', 'peak', 'noise');
+% sim_in = [t', awgn(in, noise, 'measured')];
+sim_in = [t', in];
+data = runSim('RC', sim_in);
+save('resRC_sine.mat', 'data', 'freq', 'peak', 'resolution');
 
-% combined constant + sine
-in = bias + sin(2*pi*freq*t') * peak;
+% combined ramp + sine + noise
+in = (t' * peak / duration) + sin(2*pi*freq*t') * peak;
 sim_in = [t', awgn(in, noise, 'measured')];
-datRC = runSim('RC', sim_in);
-save('resRC_comb.mat', 'datRC', 'freq', 'peak', 'bias', 'noise');
+data = runSim('RC', sim_in);
+save('resRC_comb.mat', 'data', 'freq', 'peak', 'resolution');
\ No newline at end of file
+
+clear data freq peak N t duration resolution in sim_in
\ No newline at end of file

PhysicalModels/runSim.m

@@ -11,6 +11,9 @@ stop_t = sim_in(end,1);
 steps = diff(sim_in(:, 1));
 minStep = min(steps(:));
 set_param(model, 'StopTime', num2str(stop_t), 'MaxStep', num2str(minStep));
+% set_param(model, 'StopTime', num2str(stop_t));
 sim(model);
-results = [sim_out.Time sim_out.Data];
+% results = [sim_out.Time sim_out.Data];
+[y, t] = resample(sim_out.Data, sim_out.Time, 1 / minStep);
+results = [t, y];
 end

docs/0_background.md

+---
+title: Background
+order: 0
+---
+
+System modeling can be categorized into three broad categories:
+
+## White-box models
+
+In white-box systems, the internal dynamics are evident. The system can be understood in terms of inputs, outputs, and exact internal dynamics. A white-box model assumes perfect knowledge of the system.
+
+## Grey-box models
+
+Grey-box models add some obfuscation to the system dynamics. While the structure or mathematical form of the dynamics may be known, the exact parameter values are not fullt evident. Grey-box models are constrained by the *apriori* knowledge of the system and estimate parameters within those bounds.
+
+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.
+
+## 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.
+
+
+### Further reading:
+
+* [System Identification][1]
+* [Grey box model][2]
+
+[1]: https://en.wikipedia.org/wiki/System_identification
+[2]: https://en.wikipedia.org/wiki/Grey_box_model
\ No newline at end of file

docs/1_RL_Circuit.md

+---
+title: RC Circuit Model
+order: 1
+---
+
+# RC Circuit
+
+The Resistor-Capacitor circuit represents a linear first order differential equation. Mathematically:
+
+```
+i_c = C * dV_c/dt
+
+i_r = V_r / R
+```
+
+Where
+
+* `i_c` and `i_r` are the currents flowing through the capacitor and resistor,
+* `C` is the capacitance,
+* `R` is the resistance,
+* `V_c` and `'V_r` are the voltages across the capacitor and resistor,
+* `t` is time.
+
+To make the model more realistic, a parasitic capacitance has been added in parallel to the resistor. This simulates electrical effects of conductors close by such that their electrical fields interact and affect current passing through them.
+
+![Schematic](img/1_RC_schematic.png)
+
+## Models
+
+For this model, the inputs are the driving voltage. The measurements to be simulated/learned are the current through the circuit. Known parameters are capacitance and resistance. Parasitic effects are unaccounted for and must be learned from the data.
+
+This is a first-order model. Therefore one parameter (Capacitance) is kept fixed while the other is free to be learned from the data.
+
+Model parameters are learned from input/output measurements. In this case, the input is a ramped sine wave.
+
+![Input Signal](img/1_input_signal.png)
+
+### Grey box model
+
+Using [Kirchoff's voltage law][2] and the above equations, the total voltage across the circuit can be written as:
+
+```
+V_s = 1/C ∫ i dt + i R
+
+i = [- 1 / (R C)] ∫ i dt + [1 / R] V_s
+```
+
+Using the [modern control approach][1] and taking as input (`u`), state (`x`), and output (`y`) variables:
+
+```
+x = ∫ i dt, x' = i, u = V_s
+y = x
+```
+
+We get the state-space equation:
+
+```
+x' = [- 1 / (R C)] x + [1 / R] u
+y = [- 1 / (R C)] x + [1 / R] u
+```
+
+Which yeilds the state-space matrices:
+
+```
+A = C = [- 1 / (R C)]
+B = D = [1 / R]
+```
+
+We set `R` as the learnable parameter under various conditions (noisy inputs, with/without parasitic effects). The state-space model tries to fit observed measurements with its estimate of the parameter value. The initial estimate is `R = 500`.
+
+| SNR\Parasitic 	| Yes    	| No    	|
+|---------------	|--------	|-------	|
+| 25            	| 0.0034 	| 99.99 	|
+| 250           	| 93.77  	| 99.99 	|
+
+*Table: `SNR` is the signal-to-noise ratio in dB. Parasitic effects are simulated by a 1uF capacitance in parallel to the resistor. Values are percentage of fit with the training data.*
+
+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.
+
+[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

docs/index.md

+---
+title: Documentation - Home
+order: -1
+---
+
 # Data Driven Models
 
 This project explores how combinations of physical and data based models fare against actual systems.
@@ -12,8 +17,21 @@ Data driven models do not depend on an underlying model of the system. They simp
 
 * Make a composite additive model based on physical and machine learning sub-models.
 
+## Requirements
+
+* MatLab/Simulink
+* System identification toolbox
+* Communication systems toolbox
+* SimScape library
+
 ## Project Structure
 
 ### `PhysicalModels`
 
-This directory contains files that simulate physical systems and return input/output readings. MatLab Simulink/Simscape are used for this purpose.
+This directory contains files that simulate physical systems and return input/output readings. MatLab Simulink/Simscape are used for this purpose.
\ No newline at end of file
+
+## Contents
+
+* [Background](0_background.md)
+
+* [Model: Resistor-Capacitor Circuit](1_RL_Circuit.md)
\ No newline at end of file

greyBoxRC.m

+% Estimates parameters of an RC circuit. Produces an
+% `idgrey` instance called `sys` representing the 
+% learned system parameters.
+%
+% See: https://www.mathworks.com/help/ident/ug/estimating-linear-grey-box-models.html
+%
+% Initial estimates of circuit parameters.
+% The number of parameters should match the order
+% of the system. That is, 2 parameters for a first
+% order system cannot be solved.
+par = {'resistance',500; 'capacitance', 1e-3};
+% A grey box system with identifiable parameters.
+sys_init = idgrey(@ssRC, par, 'c');
+sys_init.Structure.Parameters(2).Free = false;
+
+rcdat = load('resRC_comb.mat');
+t = rcdat.data(:,1);    % timestamps
+u = rcdat.data(:,2);    % input
+y = rcdat.data(:,3);    % output
+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));
+
+% state-space representation of RC circuit
+function [A, B, C, D] = ssRC(R, c, ~)
+    A = [-1/(R*c)];
+    B = [1/R];
+    C = A;
+    D = B;
+end
\ No newline at end of file