add code for 2D linear time-varying system

2DLinearTimeVarying/2DTimeVarying.m

@@ -0,0 +1,196 @@
+
+%% clean workspace
+clear all
+close all
+clc
+
+%% define dynamics
+epsilon = 1e-1;
+dyn = @(t,x) ([0, 1+epsilon*t; -(1+epsilon*t),0])*x;
+
+% generate data
+dt = 1e-1;
+tspan = 0:dt:10;
+x0 = [1;0];
+[tq,xq] = ode45(dyn, tspan, x0);
+% extract snapshot pairs
+xq = xq'; tq = tq';
+x = xq(:,1:end-1); y = xq(:,2:end); t = tq(2:end);
+
+% true dynamics, eigenvalues and modes
+[n, m] = size(x);
+A = zeros(2,2,m);
+evals = zeros(2,m);
+for k = 1:m
+    A(:,:,k) = [0, 1+epsilon*t(k); -(1+epsilon*t(k)),0]; % continuous time dynamics
+    evals(:,k) = eig(A(:,:,k));
+end
+
+%% visualize snapshots
+figure
+plot(tq,xq(1,:),'k-',tq,xq(2,:),'k--','MarkerSize',12,'LineWidth',1.5)
+xlabel('Time $t$','Interpreter','latex')
+ylabel('$x_1(t), x_2(t)$','Interpreter','latex')
+fl = legend('$x_1$', '$x_2$');
+set(fl,'Interpreter','latex','Box','off','FontSize',20,'Location','best');
+box on
+set(gca,'FontSize',20,'LineWidth',1)
+
+%% run and compare different DMD algorithms
+w = 10;
+
+% batch DMD
+AbatchDMD = zeros(2,2,m);
+evalsbatchDMD = zeros(2,m);
+for k = 1:w
+    AbatchDMD(:,:,k) = expm(A(:,:,k)*dt);
+end
+tic
+% start at w+1
+for k = w+1:m
+    AbatchDMD(:,:,k) = y(:,1:k)*pinv(x(:,1:k));
+end
+batchDMDelapsedTime = toc
+for k = 1:m
+    evalA = eig(AbatchDMD(:,:,k));
+    evalsbatchDMD(:,k) = log(evalA)/dt;
+end
+
+% online DMD, rho = 1
+AonlineDMD = zeros(2,2,m);
+evalsonlineDMD = zeros(2,m);
+for k = 1:w
+    AonlineDMD(:,:,k) = expm(A(:,:,k)*dt);
+end
+% initial condition
+odmd = OnlineDMD(2,1);
+odmd.initialize(x(:,1:w),y(:,1:w));
+tic
+% start at w+1
+for k = w+1:m
+    odmd.update(x(:,k),y(:,k));
+    AonlineDMD(:,:,k) = odmd.A;
+end
+onlineDMDelapsedTime = toc
+for k = 1:m
+    evalA = eig(AonlineDMD(:,:,k));
+    evalsonlineDMD(:,k) = log(evalA)/dt;
+end
+evalsonlineDMD1 = evalsonlineDMD;
+
+% online DMD, rho = 0.95
+AonlineDMD = zeros(2,2,m);
+evalsonlineDMD = zeros(2,m);
+for k = 1:w
+    AonlineDMD(:,:,k) = expm(A(:,:,k)*dt);
+end
+% initial condition
+odmd = OnlineDMD(2,0.95);
+odmd.initialize(x(:,1:w),y(:,1:w));
+tic
+% start at w+1
+for k = w+1:m
+    odmd.update(x(:,k),y(:,k));
+    AonlineDMD(:,:,k) = odmd.A;
+end
+onlineDMDelapsedTime = toc
+for k = 1:m
+    evalA = eig(AonlineDMD(:,:,k));
+    evalsonlineDMD(:,k) = log(evalA)/dt;
+end
+evalsonlineDMD2 = evalsonlineDMD;
+
+% online DMD, rho = 0.8
+AonlineDMD = zeros(2,2,m);
+evalsonlineDMD = zeros(2,m);
+for k = 1:w
+    AonlineDMD(:,:,k) = expm(A(:,:,k)*dt);
+end
+% initial condition
+odmd = OnlineDMD(2,0.8);
+odmd.initialize(x(:,1:w),y(:,1:w));
+tic
+% start at w+1
+for k = w+1:m
+    odmd.update(x(:,k),y(:,k));
+    AonlineDMD(:,:,k) = odmd.A;
+end
+onlineDMDelapsedTime = toc
+for k = 1:m
+    evalA = eig(AonlineDMD(:,:,k));
+    evalsonlineDMD(:,k) = log(evalA)/dt;
+end
+evalsonlineDMD3 = evalsonlineDMD;
+
+% streaming DMD
+evalsstreamingDMD = zeros(2,m);
+% initial
+evalsstreamingDMD(:,1) = evals(:,1);
+% streaming DMD
+max_rank = 2;
+sdmd = StreamingDMD(max_rank);
+for k = 2:m
+    sdmd = sdmd.update(x(:,k), y(:,k));
+    [~, evalA] = sdmd.compute_modes();
+    evalsstreamingDMD(:,k) = log(evalA)/dt;
+end
+
+% mini-batch DMD
+AminibatchDMD = zeros(2,2,m);
+evalsminibatchDMD = zeros(2,m);
+for k = 1:w
+    AminibatchDMD(:,:,k) = expm(A(:,:,k)*dt);
+end
+tic
+% start at w+1
+for k = w+1:m
+    AminibatchDMD(:,:,k) = y(:,k-w+1:k)*pinv(x(:,k-w+1:k));
+end
+minibatchDMDelapsedTime = toc
+for k = 1:m
+    evalA = eig(AminibatchDMD(:,:,k));
+    evalsminibatchDMD(:,k) = log(evalA)/dt;
+end
+
+% window DMD
+AwindowDMD = zeros(2,2,m);
+evalswindowDMD = zeros(2,m);
+for k = 1:w
+    AwindowDMD(:,:,k) = expm(A(:,:,k)*dt);
+end
+% initialization
+wdmd = WindowDMD(2,w,1);
+wdmd.initialize(x(:,1:w), y(:,1:w));
+tic
+% start at w+1
+for k = w+1:m
+    wdmd.update(x(:,k), y(:,k));
+    AwindowDMD(:,:,k) = wdmd.A;
+end
+windowDMDelapsedTime = toc
+for k = 1:m
+    evalA = eig(AwindowDMD(:,:,k));
+    evalswindowDMD(:,k) = log(evalA)/dt;
+end
+
+%% Compare various DMD algorithm
+index = w+1:length(t);
+
+% plot Imaginary part
+
+figure, hold on
+plot(t,imag(evals(1,:)),'k-','LineWidth',1)
+plot(t(1,index),imag(evalsbatchDMD(1,index)),'s-','LineWidth',1,'MarkerSize',12,'MarkerIndices',1:10:length(index))
+plot(t(1,index),imag(evalsminibatchDMD(1,index)),'o-','LineWidth',1,'MarkerSize',12,'MarkerIndices',1:10:length(index))
+plot(t(1,index),imag(evalsstreamingDMD(1,index)),'x-','LineWidth',1,'MarkerSize',12,'MarkerIndices',7:10:length(index))
+plot(t(1,index),imag(evalsonlineDMD1(1,index)),'>-','LineWidth',1,'MarkerSize',12,'MarkerIndices',4:10:length(index))
+plot(t(1,index),imag(evalsonlineDMD2(1,index)),'*-','LineWidth',1,'MarkerSize',12,'MarkerIndices',1:10:length(index))
+plot(t(1,index),imag(evalsonlineDMD3(1,index)),'d-','LineWidth',1,'MarkerSize',12,'MarkerIndices',4:10:length(index))
+plot(t(1,index),imag(evalswindowDMD(1,index)),'+-','LineWidth',1,'MarkerSize',12,'MarkerIndices',7:10:length(index))
+xlabel('Time $t$','Interpreter','latex'), ylabel('Im$(\lambda)$','Interpreter','latex')
+fl = legend('true','batch','mini-batch','streaming','online, $\rho=1$','online, $\rho=0.95$','online, $\rho=0.8$','window');
+set(fl,'Interpreter','latex','Location','northwest','FontSize',20,'Box','off');
+xlim([0,10]), ylim([1,2])
+
+box on
+set(gca,'FontSize',20,'LineWidth',1)

2DLinearTimeVarying/OnlineDMD.m

@@ -0,0 +1,129 @@
+% OnlineDMD is a class that implements online dynamic mode decomposition
+% The time complexity (multiply-add operation for one iteration) is O(4n^2)
+% , and space complexity is O(2n^2), where n is the state dimension.
+%
+% Algorithm description:
+%       At time step k, define two matrix X(k) = [x(1),x(2),...,x(k)], 
+%       Y(k) = [y(1),y(2),...,y(k)], that contain all the past snapshot 
+%       pairs, where x(k), y(k) are the n dimensional state vector, 
+%       y(k) = f(x(k)) is the image of x(k), f() is the dynamics. 
+%       Here, if the (discrete-time) dynamics are given by z(k) = f(z(k-1))
+%       , then x(k), y(k) should be measurements corresponding to 
+%       consecutive states z(k-1) and z(k).
+%       At time step k+1, we need to include new snapshot pair x(k+1), y(k+1)
+%       We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively 
+%       by efficient rank-1 updating online DMD algorithm.
+%       An exponential weighting factor can be used to place more weight on
+%       recent data.
+%
+% Usage:
+%       odmd = OnlineDMD(n,weighting)
+%       odmd.initialize(Xq,Yq)
+%       odmd.initilizeghost()
+%       odmd.update(x,y)
+%       [evals, modes] = odmd.computemodes()
+%        
+% properties:
+%       n: state dimension
+%       weighting: weighting factor in (0,1]
+%       timestep: number of snapshot pairs processed
+%       A: DMD matrix, size n by n
+%       P: matrix that contains information about past snapshots, size n by n
+% 
+% methods:
+%       initialize(Xq, Yq), initialize online DMD algorithm with q snapshot
+%                           pairs stored in (Xq, Yq)
+%       initializeghost(), initialize online DMD algorithm with epsilon 
+%                           small (1e-15) ghost snapshot pairs before t=0
+%       update(x,y), update when new snapshot pair (x,y) becomes available
+%                   Here, if the (discrete-time) dynamics are given by 
+%                   z(k) = f(z(k-1)), then (x,y) should be measurements 
+%                   correponding to consecutive states z(k-1) and z(k).
+%       computemodes(), compute and return DMD eigenvalues and DMD modes
+%
+% Authors: 
+%   Hao Zhang
+%   Clarence W. Rowley
+% 
+% Reference:
+% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta,
+% ``Online Dynamic Mode Decomposition for Time-varying Systems", 
+% in production, 2017. Available on arXiv.
+% 
+% Created:
+%   April 2017.
+%
+% To look up the documentation in the command window, type help OnlineDMD
+
+classdef OnlineDMD < handle
+    properties
+        n = 0;                      % state dimension
+        weighting = 1;                 % weighting factor
+        timestep = 0;               % number of snapshots processed
+        A;          % DMD matrix
+        P;          % matrix that contains information about past snapshots
+    end
+
+    methods
+        function obj = OnlineDMD(n,weighting)
+            % Creat an object for online DMD
+            % Usage: odmd = OnlineDMD(n,weighting)
+            if nargin == 2
+                obj.n = n;
+                obj.weighting = weighting;
+                obj.A = zeros(n,n);
+                obj.P = zeros(n,n);
+            end
+        end
+
+        function initialize(obj, Xq, Yq)
+            % Initialize OnlineDMD with q snapshot pairs stored in (Xq, Yq)
+            % Usage: odmd.initialize(Xq,Yq)
+            q = length(Xq(1,:));
+            if(obj.timestep == 0 && rank(Xq) == obj.n)
+                weight = (sqrt(obj.weighting)).^(q-1:-1:0);
+                Xq = Xq.*weight;
+                Yq = Yq.*weight;
+                obj.A = Yq*pinv(Xq);
+                obj.P = inv(Xq*Xq')/obj.weighting;
+            end
+            obj.timestep = obj.timestep + q;
+        end
+
+        function initializeghost(obj)
+            % Initialize online DMD with epsilon small (1e-15) ghost 
+            % snapshot pairs before t=0
+            % Usage: odmd.initilizeghost()
+            epsilon = 1e-15;
+            obj.A = randn(obj.n, obj.n);
+            obj.P = (1/epsilon)*eye(obj.n);
+        end
+
+        function update(obj, x, y)
+            % Update the DMD computation with a new pair of snapshots (x,y)
+            % Here, if the (discrete-time) dynamics are given by z(k) = 
+            % f(z(k-1)), then (x,y) should be measurements correponding to 
+            % consecutive states z(k-1) and z(k).
+            % Usage: odmd.update(x, y)
+
+            % compute P*x matrix vector product beforehand
+            Px = obj.P*x;
+            % Compute gamma
+            gamma = 1/(1+x'*Px);
+            % Update A
+            obj.A = obj.A + (gamma*(y-obj.A*x))*Px';
+            % Update P, group Px*Px' to ensure positive definite
+            obj.P = (obj.P - gamma*(Px*Px'))/obj.weighting;
+            % ensure P is SPD by taking its symmetric part
+            obj.P = (obj.P+(obj.P)')/2;
+            % time step + 1
+            obj.timestep = obj.timestep + 1;
+        end
+
+        function [evals, modes] = computemodes(obj)
+            % Compute DMD eigenvalues and DMD modes at current time
+            % Usage: [evals, modes] = odmd.modes()
+            [modes, evals] = eig(obj.A, 'vector');
+        end
+    end
+end