diff --git a/README.md b/README.md index e78e3fb..f0d1fa7 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,15 @@ # dmdtools A library of tools for computing variants of Dynamic Mode Decomposition + +## algorithms implemented in Matlab +1. Standard batch processed total-least-squares DMD (TDMD) +2. Streaming DMD (including total-least-squares) +3. Online DMD (including weighted version) +4. Window DMD (including weighted version) + +## algorithms implemented in Python +1. Standard batch processed DMD +2. Kernel DMD with polynomial kernel +3. Streaming DMD +4. Online DMD (including weighted version) +5. Window DMD (including weighted version) \ No newline at end of file diff --git a/matlab/OnlineDMD.m b/matlab/OnlineDMD.m new file mode 100644 index 0000000..9d82f6e --- /dev/null +++ b/matlab/OnlineDMD.m @@ -0,0 +1,130 @@ +% 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," +% arXiv preprint arXiv:1707.02876, 2017. +% +% 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 in (0,1] + 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.timestep = 0; + 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 \ No newline at end of file diff --git a/matlab/README b/matlab/README deleted file mode 100755 index c448a34..0000000 --- a/matlab/README +++ /dev/null @@ -1,32 +0,0 @@ -README for a Streaming Dynamic Mode Decomposition (SDMD) -======================================================== - -The Matlab scripts contained in this distribution contain: (1) a Matlab class for an -object-oriented impelentation of incrementally updated DMD, and (2) an example -to demonstrate the use of this class on an arbitrary dynamical system. These scripts -are supplementary materials associated with Hemati et al. (2014). - - -SDMD Distribution Contents -========================== - -StreamingDMD.m --------------- -A Matlab class that implements SDMD in an object-oriented fashion, -based on the formulation in Hemati et al. (2014). - -sdmd_run.m ----------- -Demonstrates usage of the StreamingDMD class in solving an arbitrary -dynamical system with two characteristic frequencies. - - -Note: use the `help' and `doc' commands within Matlab to access more detailed -descriptions of these scripts. - - -References -========== -Maziar S. Hemati, Matthew O. Williams, and Clarence W. Rowley, -``Dynamic Mode Decomposition for Large and Streaming Datasets,'' -Physics of Fluids (2014). \ No newline at end of file diff --git a/matlab/README.md b/matlab/README.md new file mode 100644 index 0000000..820dfea --- /dev/null +++ b/matlab/README.md @@ -0,0 +1,15 @@ +# README for Matlab implementation of variants of Dynamic Mode Decomposition + +## Implementation +1.**tdmd.m** implements the total-least-squares DMD function. +2.**StreamingDMD.m** implements **StramingDMD** class. +3.**StreamingTDMD.m** implements **StreamingTDMD** class. +4.**OnlineDMD.m** implements **OnlineDMD** class. +5.**WindomDMD.m** implements **WindowDMD** class. + +## Demos +1.**tdmd_run.m** demostrates the use of total-least-squares DMD function. +2.**sdmd_run.m** demostrates the use of **StreamingDMD** class. +3.**stdmd_run.m** demostrates the use of **StreamingTDMD** class. +4.**online_demo.m** demostrates the use of **OnlineDMD** class. +5.**window_demo.m** demostrates the use of **WindowDMD** class. \ No newline at end of file diff --git a/matlab/WindowDMD.m b/matlab/WindowDMD.m new file mode 100644 index 0000000..75b95a3 --- /dev/null +++ b/matlab/WindowDMD.m @@ -0,0 +1,156 @@ +% WindowDMD is a class that implements window dynamic mode decomposition +% The time complexity (multiply-add operation for one iteration) is +% O(8n^2), and space complexity is O(2wn+2n^2), where n is the state +% dimension, w is the window size. +% +% Algorithm description: +% At time step k, define two matrix +% X(k) = [x(k-w+1),x(k-w+2),...,x(k)], Y(k) = [y(k-w+1),y(k-w+2), +% ...,y(k)], that contain the recent w snapshot pairs from a finite +% time window, 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 forget old snapshot pair xold = +% x(k-w+1), yold = y(k-w+1), and remember new snapshot pair xnew = +% x(k+1), ynew = y(k+1). +% We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively +% by efficient rank-2 updating window DMD algorithm. +% An exponential weighting factor can be used to place more weight on +% recent data. +% +% Usage: +% wdmd = WindowDMD(n,w,weighting) +% wdmd.initialize(Xw,Yw) +% wdmd.update(xnew, ynew) +% [evals, modes] = wdmd.computemodes() +% +% properties: +% n: state dimension +% w: finite time window size +% weighting: weighting factor in (0,1] +% timestep: number of snapshot pairs processed +% Xw: recent w snapshots x stored in Xw, size n by w +% Yw: recent w snapshots y stored in Yw, size n by w +% A: DMD matrix for w snapshot pairs, size n by n +% P: Matrix that contains information about recent w snapshots, +% size n by n +% +% methods: +% initialize(Xw, Yw), initialize window DMD algorithm +% update(xnew, ynew), update by forgetting old snapshot pairs, +% and remeber new snapshot pair, move sliding window forward +% At time k+1, X(k+1) = [x(k-w+2),x(k-w+2),...,x(k+1)], +% Y(k+1) = [y(k-w+2),y(k-w+2),...,y(k+1)], +% we should take xnew = x(k+1), ynew = y(k+1) +% 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," +% arXiv preprint arXiv:1707.02876, 2017. +% +% Created: +% April 2017. +% +% To look up the documentation, type help WindowDMD + +classdef WindowDMD < handle + properties + n = 0; % state dimension + w = 0; % window size + weighting = 1; % weighting factor in (0,1] + timestep = 0; % number of snapshots processed + Xw; % recent w snapshots x stored in matrix Xw + Yw; % recent w snapshots y stored in matrix Yw + A; % DMD matrix for w snapshot pairs, size n by n + P; % Matrix that contains information about recent w snapshots + % , size n by n + end + + methods + function obj = WindowDMD(n, w, weighting) + % Creat an object for window DMD + % Usage: wdmd = WindowDMD(n,w,weighting) + if nargin == 3 + obj.n = n; + obj.w = w; + obj.weighting = weighting; + obj.timestep = 0; + obj.Xw = zeros(n,w); + obj.Yw = zeros(n,w); + obj.A = zeros(n,n); + obj.P = zeros(n,n); + end + end + + function initialize(obj, Xw, Yw) + % Initialize WindowDMD with w snapshot pairs stored in (Xw, Yw) + % Usage: wdmd.initialize(Xw,Yw) + + % initialize Xw, Yw + obj.Xw = Xw; obj.Yw = Yw; + % initialize A, P + q = length(Xw(1,:)); + if(obj.timestep == 0 && obj.w == q && rank(Xw) == obj.n) + weight = (sqrt(obj.weighting)).^(q-1:-1:0); + Xw = Xw.*weight; + Yw = Yw.*weight; + obj.A = Yw*pinv(Xw); + obj.P = inv(Xw*Xw')/obj.weighting; + end + obj.timestep = obj.timestep + q; + end + + function update(obj, xnew, ynew) + % Update the DMD computation by sliding the finite time window + % forward. + % Forget the oldest pair of snapshots (xold, yold), and + % remembers the newest pair of snapshots (xnew, ynew) in the + % new time window. If the new finite time window at time step + % k+1 includes recent w snapshot pairs as + % X(k+1) = [x(k-w+2),x(k-w+3),...,x(k+1)], + % Y(k+1) = [y(k-w+2),y(k-w+3),...,y(k+1)], + % where y(k) = f(x(k)) and f is the dynamics, then we should + % take xnew = x(k+1), ynew = y(k+1) + % Usage: wdmd.update(xnew, ynew) + + % define old snapshots to be discarded + xold = obj.Xw(:,1); yold = obj.Yw(:,1); + % Update recent w snapshots + obj.Xw = [obj.Xw(:,2:end), xnew]; + obj.Yw = [obj.Yw(:,2:end), ynew]; + + % direct rank-2 update + % define matrices + U = [xold, xnew]; V = [yold, ynew]; + C = diag([-(obj.weighting)^(obj.w),1]); + % compute PkU matrix matrix product beforehand + PkU = obj.P*U; + % compute AkU matrix matrix product beforehand + AkU = obj.A*U; + % compute Gamma + Gamma = inv(inv(C)+U'*PkU); + % update A + obj.A = obj.A + (V-AkU)*(Gamma*PkU'); + % update P + obj.P = (obj.P - PkU*(Gamma*PkU'))/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 step + % Usage: [evals, modes] = wdmd.computemodes() + [modes, evals] = eig(obj.A, 'vector'); + end + end +end \ No newline at end of file diff --git a/matlab/online_demo.m b/matlab/online_demo.m new file mode 100644 index 0000000..d2fa825 --- /dev/null +++ b/matlab/online_demo.m @@ -0,0 +1,117 @@ +% An example to demonstrate online dynamic mode decomposition +% +% We take a 2D time varying system given by dx/dt = A(t)x +% where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +% w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +% are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. +% +% 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, +% we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force +% batch DMD, and by efficient rank-1 updating online DMD algrithm. Batch DMD +% computes DMD matrix by brute-force taking the pseudo-inverse directly. +% Online DMD computes the DMD matrix by using efficient rank-1 update idea. +% +% We compare the performance of online DMD (with weighting=1,0.9) with the +% brute-force batch DMD approach in terms of tracking time varying eigenvalues, +% by comparison with the analytical solution. Online DMD (weighting=1) and +% batch DMD should agree with each other (up to machine round-offer errors). +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% References: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% "Online Dynamic Mode Decomposition for Time-varying Systems," +% arXiv preprint arXiv:1707.02876, 2017. +% +% Date created: April 2017 + +% 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); time = tq(2:end); +% true dynamics, eigenvalues +[n, m] = size(x); +A = zeros(n,n,m); +evals = zeros(n,m); +for k = 1:m + A(:,:,k) = [0, 1+epsilon*time(k); -(1+epsilon*time(k)),0]; % continuous time dynamics + evals(:,k) = eig(A(:,:,k)); % analytical continuous time eigenvalues +end + + +% visualize snapshots +figure, hold on +plot(tq,xq(1,:),'x-',tq,xq(2,:),'o-','LineWidth',2) +xlabel('Time','Interpreter','latex') +title('Snapshots','Interpreter','latex') +fl = legend('$x_1(t)$','$x_2(t)$'); +set(fl,'Interpreter','latex'); +box on +set(gca,'FontSize',20,'LineWidth',2) + + +% batch DMD +q = 10; +AbatchDMD = zeros(n,n,m); +evalsbatchDMD = zeros(n,m); +tic +for k = q+1:m + AbatchDMD(:,:,k) = y(:,1:k)*pinv(x(:,1:k)); + evalsbatchDMD(:,k) = log(eig(AbatchDMD(:,:,k)))/dt; +end +elapsed_time = toc; +fprintf('Batch DMD, elapsed time: %f seconds\n', elapsed_time) + +% Online DMD weighting = 1 +evalsonlineDMD1 = zeros(n,m); +% creat object and initialize with first q snapshot pairs +odmd = OnlineDMD(n,1); +odmd.initialize(x(:,1:q),y(:,1:q)); +% online DMD +tic +for k = q+1:m + odmd.update(x(:,k),y(:,k)); + evalsonlineDMD1(:,k) = log(eig(odmd.A))/dt; +end +elapsed_time = toc; +fprintf('Online DMD, weighting = 1, elapsed time: %f seconds\n', elapsed_time) + +% Online DMD, weighting = 0.9 +evalsonlineDMD2 = zeros(n,m); +% creat object and initialize with first q snapshot pairs +odmd = OnlineDMD(n,0.9); +odmd.initialize(x(:,1:q),y(:,1:q)); +% online DMD +tic +for k = q+1:m + odmd.update(x(:,k),y(:,k)); + evalsonlineDMD2(:,k) = log(eig(odmd.A))/dt; +end +elapsed_time = toc; +fprintf('Online DMD, weighting = 0.9, elapsed time: %f seconds\n', elapsed_time) + + +% visualize imaginary part of the continous time eigenvalues +% from true, batch, online (rho=1), and online (rho=0.9) +updateindex = q+1:m; +figure, hold on +plot(time,imag(evals(1,:)),'k-','LineWidth',2) +plot(time(updateindex),imag(evalsbatchDMD(1,updateindex)),'-','LineWidth',2) +plot(time(updateindex),imag(evalsonlineDMD1(1,updateindex)),'--','LineWidth',2) +plot(time(updateindex),imag(evalsonlineDMD2(1,updateindex)),'--','LineWidth',2) +xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') +fl = legend('True','Batch','Online, $wf=1$','Online, $wf=0.9$'); +set(fl,'Interpreter','latex','Location','northwest'); +ylim([1,2]), xlim([0,10]) +box on +set(gca,'FontSize',20,'LineWidth',2) \ No newline at end of file diff --git a/matlab/window_demo.m b/matlab/window_demo.m new file mode 100644 index 0000000..2f81830 --- /dev/null +++ b/matlab/window_demo.m @@ -0,0 +1,121 @@ +% An example to demonstrate window dynamic mode decomposition +% +% We take a 2D time varying system given by dx/dt = A(t)x +% where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +% w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +% are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. +% +% At time step k, define two matrix X(k) = [x(k-w+1),x(k-w+2),...,x(k)], +% Y(k) = [y(k-w+1),y(k-w+2),...,y(k)], that contain the recent w snapshot pairs +% from a finite time window, we would like to compute Ak = Yk*pinv(Xk). This can +% be done by brute-force mini-batch DMD, and by efficient rank-2 updating window +% DMD algrithm. For window DMD, at time k+1, we need to forget the old snapshot +% pair xold = x(k-w+1), yold = y(k-w+1), and remember the new snapshot pair xnew +% = x(k+1), ynew = y(k+1). Mini-batch DMD computes DMD matrix by taking the +% pseudo-inverse directly. Window DMD computes the DMD matrix by using efficient +% rank-2 update idea. +% +% We compare the performance of window DMD with the brute-force mini-batch DMD +% approach in terms of tracking time varying eigenvalues, by comparison with +% the analytical solution. They should agree with each other (up to machine +% round-offer errors). +% +% Authors: +% Hao Zhang +% Clarence W. Rowley +% +% References: +% Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, +% "Online Dynamic Mode Decomposition for Time-varying Systems," +% arXiv preprint arXiv:1707.02876, 2017. +% +% Date created: April 2017 + +% 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); time = tq(2:end); +% true dynamics, eigenvalues +[n, m] = size(x); +A = zeros(n,n,m); +evals = zeros(n,m); +for k = 1:m + A(:,:,k) = [0, 1+epsilon*time(k); -(1+epsilon*time(k)),0]; % continuous time dynamics + evals(:,k) = eig(A(:,:,k)); % analytical continuous time eigenvalues +end + + +% visualize snapshots +figure, hold on +plot(tq,xq(1,:),'x-',tq,xq(2,:),'o-','LineWidth',2) +xlabel('Time','Interpreter','latex') +title('Snapshots','Interpreter','latex') +fl = legend('$x_1(t)$','$x_2(t)$'); +set(fl,'Interpreter','latex'); +box on +set(gca,'FontSize',20,'LineWidth',2) + + +% mini-batch DMD +w = 10; % storage time window size, store recent w snapshot pairs +AminibatchDMD = zeros(n,n,m); +evalsminibatchDMD = zeros(n,m); +% mini-batch DMD +tic +for k = w+1:m + AminibatchDMD(:,:,k) = y(:,k-w+1:k)*pinv(x(:,k-w+1:k)); + evalsminibatchDMD(:,k) = log(eig(AminibatchDMD(:,:,k)))/dt; +end +elapsed_time = toc; +fprintf('Mini-batch DMD, w=10, elapsed time: %f seconds\n', elapsed_time) + + +% window DMD, weighting = 1 +evalswindowDMD1 = zeros(n,m); +% creat object and initialize with first w snapshot pairs +wdmd = WindowDMD(n,w,1); +wdmd.initialize(x(:,1:w), y(:,1:w)); +% window DMD +tic +for k = w+1:m + wdmd.update(x(:,k), y(:,k)); + evalswindowDMD1(:,k) = log(eig(wdmd.A))/dt; +end +elapsed_time = toc; +fprintf('Window DMD, w = 10, weighting = 1, elapsed time: %f seconds\n', elapsed_time) + +% window DMD, weighting = 0.5 +evalswindowDMD2 = zeros(n,m); +% creat object and initialize with first w snapshot pairs +wdmd = WindowDMD(n,w,0.5); +wdmd.initialize(x(:,1:w), y(:,1:w)); +% window DMD +tic +for k = w+1:m + wdmd.update(x(:,k), y(:,k)); + evalswindowDMD2(:,k) = log(eig(wdmd.A))/dt; +end +elapsed_time = toc; +fprintf('Window DMD, w = 10, weighting = 0.5, elapsed time: %f seconds\n', elapsed_time) + +% visualize imaginary part of the continous time eigenvalues +% from true, mini-batch, and window +updateindex = w+1:m; +figure, hold on +plot(time,imag(evals(1,:)),'k-','LineWidth',2) +plot(time(updateindex),imag(evalsminibatchDMD(1,updateindex)),'-','LineWidth',2) +plot(time(updateindex),imag(evalswindowDMD1(1,updateindex)),'--','LineWidth',2) +plot(time(updateindex),imag(evalswindowDMD2(1,updateindex)),'--','LineWidth',2) +xlabel('Time','Interpreter','latex'), ylabel('Im($\lambda_{DMD}$)','Interpreter','latex') +fl = legend('True','Mini-batch, $w=10$','Window, $w=10$, $wf=1$','Window, $w=10$, $wf=0.5$'); +set(fl,'Interpreter','latex','Location','northwest'); +ylim([1,2]), xlim([0,10]) +box on +set(gca,'FontSize',20,'LineWidth',2) \ No newline at end of file diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000..ac4f9a6 --- /dev/null +++ b/python/README.md @@ -0,0 +1,18 @@ +# README for Python implementation of variants of Dynamic Mode Decomposition + +## folders +1. **dmdtools** folder contains implementation of variants of Dynamic Mode Decomposition +2. **scripts** folder contains demos of various DMD algrithm +3. **tests** folder contains tests for DMD class and kernel DMD class. + +## implementations contained in **dmdtools** folder +1. **batch.py** implements Standard batch processed DMD (including total-least-squares) and Kernel DMD with polynomial kernel. +2. **streaming.py** implements Streaming DMD. +3. **online.py** implements online DMD (including weighted online DMD). +4. **window.py** implements window DMD. + +## demostrations contained in **scripts** folder +1. **total_dmd_example.py** demostrates the use of DMD and total-least-squares DMD. +2. **streaming_dmd_example.py** demostrates the use of Streaming DMD. +3. **online_demo.py** demostrates the use of online DMD +4. **window_demo.py** demostrates the use of window DMD \ No newline at end of file diff --git a/python/dmdtools/__init__.py b/python/dmdtools/__init__.py index b6635bc..4aac6ee 100644 --- a/python/dmdtools/__init__.py +++ b/python/dmdtools/__init__.py @@ -1,4 +1,4 @@ from batch import * -#from regularized import * from streaming import * - +from online import * +from window import * \ No newline at end of file diff --git a/python/dmdtools/online.py b/python/dmdtools/online.py new file mode 100644 index 0000000..fd79138 --- /dev/null +++ b/python/dmdtools/online.py @@ -0,0 +1,120 @@ +# -*- coding: utf-8 -*- +import numpy as np + + +class OnlineDMD: + """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 correponding to consecutive + states z(k-1) and z(k). + We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively + by efficient rank-1 updating online DMD algrithm. + 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 (i.e., current time step) + 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 first 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 DMD computation when new snapshot pair (x,y) + becomes available + computemodes(), compute and return DMD eigenvalues and DMD modes + + Authors: + Hao Zhang + Clarence W. Rowley + + References: + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + "Online Dynamic Mode Decomposition for Time-varying Systems," + arXiv preprint arXiv:1707.02876, 2017. + + Date created: April 2017 + + To import the OnlineDMD class, add import online at head of Python scripts. + To look up this documentation, type help(online.OnlineDMD) or + online.OnlineDMD? + """ + def __init__(self, n=0, weighting=1): + """ + Creat an object for online DMD + Usage: odmd = OnlineDMD(n,weighting) + """ + self.n = n + self.weighting = weighting + self.timestep = 0 + self.A = np.zeros([n,n]) + self.P = np.zeros([n,n]) + + def initialize(self, Xq, Yq): + """Initialize online DMD with first q snapshot pairs stored in (Xq, Yq) + Usage: odmd.initialize(Xq,Yq) + """ + q = len(Xq[0,:]) + Xqhat, Yqhat = np.zeros(Xq.shape), np.zeros(Yq.shape) + if self.timestep == 0 and np.linalg.matrix_rank(Xq) == self.n: + weight = np.sqrt(self.weighting)**range(q-1,-1,-1) + Xqhat, Yqhat = weight*Xq, weight*Yq + self.A = Yqhat.dot(np.linalg.pinv(Xqhat)) + self.P = np.linalg.inv(Xqhat.dot(Xqhat.T))/self.weighting + self.timestep += q + + def initializeghost(self): + """Initialize online DMD with epsilon small (1e-15) ghost snapshot pairs + before t=0 + Usage: odmd.initilizeghost() + """ + epsilon=1e-15 + alpha = 1.0/epsilon + self.A = np.random.randn(self.n, self.n) + self.P = alpha*np.identity(self.n) + + def update(self, 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 = self.P.dot(x) + # compute gamma + gamma = 1.0/(1 + x.T.dot(Px)) + # update A + self.A += np.outer(gamma*(y-self.A.dot(x)),Px) + # update P, group Px*Px' to ensure positive definite + self.P = (self.P - gamma*np.outer(Px,Px))/self.weighting + # ensure P is SPD by taking its symmetric part + self.P = (self.P + self.P.T)/2 + # time step + 1 + self.timestep += 1 + + def computemodes(self): + """Compute and return DMD eigenvalues and DMD modes at current time step + Usage: evals, modes = odmd.computemodes() + """ + evals, modes = np.linalg.eig(self.A) + return evals, modes \ No newline at end of file diff --git a/python/dmdtools/window.py b/python/dmdtools/window.py new file mode 100644 index 0000000..2d072fd --- /dev/null +++ b/python/dmdtools/window.py @@ -0,0 +1,134 @@ +# -*- coding: utf-8 -*- +import numpy as np + + +class WindowDMD: + """WindowDMD is a class that implements window dynamic mode decomposition + The time complexity (multiply–add operation for one iteration) is O(8n^2), + and space complexity is O(2wn+2n^2), where n is the state dimension, w is + the window size. + + Algorithm description: + At time step k, define two matrix X(k) = [x(k-w+1),x(k-w+2),...,x(k)], + Y(k) = [y(k-w+1),y(k-w+2),...,y(k)], that contain the recent w snapshot + pairs from a finite time window, 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 correponding to consecutive + states z(k-1) and z(k). + At time k+1, we need to forget the old snapshot pair xold = x(k-w+1), + yold = y(k-w+1), and remember the new snapshot pair xnew = x(k+1), + ynew = y(k+1). + We would like to update the DMD matrix Ak = Yk*pinv(Xk) recursively + by efficient rank-2 updating window DMD algrithm. + An exponential weighting factor can be used to place more weight on + recent data. + + Usage: + wdmd = WindowDMD(n,w) + wdmd.initialize(Xw,Yw) + wdmd.update(xnew,ynew) + evals, modes = wdmd.computemodes() + + properties: + n: state dimension + w: window size + weighting: weighting factor in (0,1] + timestep: number of snapshot pairs processed (i.e., current time step) + Xw: recent w snapshots x stored in Xw, size n by w + Yw: recent w snapshots y stored in Yw, size n by w + A: DMD matrix, size n by n + P: Matrix that contains information about recent w snapshots, size n by n + + methods: + initialize(Xw, Yw), initialize window DMD algorithm with w snapshot pairs + update(xnew,ynew), update DMD computation by adding a new snapshot pair + computemodes(), compute and return DMD eigenvalues and DMD modes + + Authors: + Hao Zhang + Clarence W. Rowley + + References: + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + "Online Dynamic Mode Decomposition for Time-varying Systems," + arXiv preprint arXiv:1707.02876, 2017. + + Date created: April 2017 + + To import the WindowDMD class, add import window at head of Python scripts. + To look up this documentation, type help(window.WindowDMD) + or window.WindowDMD? + """ + def __init__(self, n=0, w=0, weighting=1): + """ + Creat an object for window DMD + Usage: wdmd = WindowDMD(n,w,weighting) + """ + self.n = n + self.w = w + self.weighting = weighting + self.timestep = 0 + self.Xw = np.zeros([n,w]) + self.Yw = np.zeros([n,w]) + self.A = np.zeros([n,n]) + self.P = np.zeros([n,n]) + + def initialize(self, Xw, Yw): + """Initialize window DMD with first w snapshot pairs stored in (Xw, Yw) + Usage: wdmd.initialize(Xw,Yw) + """ + # initialize Xw, Yw + self.Xw, self.Yw = Xw, Yw + # initialize A, P + q = len(Xw[0,:]) + if self.timestep == 0 and self.w == q \ + and np.linalg.matrix_rank(Xw) == self.n: + weight = np.sqrt(self.weighting)**range(q-1,-1,-1) + Xwhat, Ywhat = weight*Xw, weight*Yw + self.A = Ywhat.dot(np.linalg.pinv(Xwhat)) + self.P = np.linalg.inv(Xwhat.dot(Xwhat.T))/self.weighting + self.timestep += q + + def update(self, xnew, ynew): + """Update the DMD computation by sliding the finite time window forward + Forget the oldest pair of snapshots (xold, yold), and remembers the newest + pair of snapshots (xnew, ynew) in the new time window. If the new finite + time window at time step k+1 includes recent w snapshot pairs as + X(k+1) = [x(k-w+2),x(k-w+3),...,x(k+1)], Y(k+1) = [y(k-w+2),y(k-w+3), + ...,y(k+1)], where y(k) = f(x(k)) and f is the dynamics, then we should + take xnew = x(k+1), ynew = y(k+1) + Usage: wdmd.update(xnew, ynew) + """ + # define old snapshots to be discarded + xold, yold = self.Xw[:,0], self.Yw[:,0] + # Update recent w snapshots + self.Xw = np.column_stack((self.Xw[:,1:], xnew)) + self.Yw = np.column_stack((self.Yw[:,1:], ynew)) + + # direct rank-2 update + # define matrices + U, V = np.column_stack((xold, xnew)), np.column_stack((yold, ynew)) + C = np.diag([-(self.weighting)**(self.w),1]) + # compute PkU matrix matrix product beforehand + PkU = self.P.dot(U) + # compute AkU matrix matrix product beforehand + AkU = self.A.dot(U) + # compute Gamma + Gamma = np.linalg.inv(np.linalg.inv(C)+U.T.dot(PkU)) + # update A + self.A += (V-AkU).dot(Gamma).dot(PkU.T) + # update P + self.P = (self.P - PkU.dot(Gamma).dot(PkU.T))/self.weighting + # ensure P is SPD by taking its symmetric part + self.P = (self.P + self.P.T)/2 + + # time step + 1 + self.timestep += 1 + + def computemodes(self): + """Compute and return DMD eigenvalues and DMD modes at current time step + Usage: evals, modes = wdmd.computemodes() + """ + evals, modes = np.linalg.eig(self.A) + return evals, modes \ No newline at end of file diff --git a/python/scripts/online_demo.py b/python/scripts/online_demo.py new file mode 100644 index 0000000..6a824c4 --- /dev/null +++ b/python/scripts/online_demo.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +""" +An example to demonstrate online dynamic mode decomposition + +We take a 2D time varying system given by dx/dt = A(t)x +where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. + +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, +we would like to compute Ak = Yk*pinv(Xk). This can be done by brute-force +batch DMD, and by efficient rank-1 updating online DMD algrithm. Batch DMD +computes DMD matrix by brute-force taking the pseudo-inverse directly. +Online DMD computes the DMD matrix by using efficient rank-1 update idea. + +We compare the performance of online DMD (with weighting=1,0.9) with the +brute-force batch DMD approach in terms of tracking time varying eigenvalues, +by comparison with the analytical solution. Online DMD (weighting=1) and +batch DMD should agree with each other (up to machine round-offer errors). + +Authors: + Hao Zhang + Clarence W. Rowley + +References: + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + "Online Dynamic Mode Decomposition for Time-varying Systems," + arXiv preprint arXiv:1707.02876, 2017. + +Date created: April 2017 +""" + + +import sys +sys.path.append('..') + +import dmdtools +import numpy as np +from scipy.integrate import odeint +import time +import matplotlib.pyplot as plt + + +# define dynamics +epsilon = 1e-1 +def dyn(x,t): + x1, x2 = x + dxdt = [(1+epsilon*t)*x2,-(1+epsilon*t)*x1] + return dxdt +# integrate from initial condition [1,0] +tspan = np.linspace(0,10,101) +dt = 0.1 +x0 = [1,0] +xsol = odeint(dyn,x0,tspan).T +# extract snapshots +x, y = xsol[:,:-1], xsol[:,1:] +t = tspan[1:] +# true dynamics, true eigenvalues +n, m = len(x[:,0]), len(x[0,:]) +A = np.empty((n,n,m)) +evals = np.empty((n,m),dtype=complex) +for k in range(m): + A[:,:,k] = np.array([[0,(1+epsilon*t[k])],[-(1+epsilon*t[k]),0]]) + evals[:,k] = np.linalg.eigvals(A[:,:,k]) + + +# visualize snapshots +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(tspan, xsol[0,:], 'bs-', linewidth=2.0, label='$x_1(t)$') +plt.plot(tspan, xsol[1,:], 'g^-', linewidth=2.0, label='$x_2(t)$') +plt.legend(loc='best',fontsize=20 ,shadow=True) +plt.xlabel('Time', fontsize=20) +plt.title('Snapshots', fontsize=20) +plt.tick_params(labelsize=20) +plt.grid() +plt.show() + + +# batch DMD +q = 10 +AbatchDMD = np.empty((n,n,m)) +evalsbatchDMD = np.empty((n,m),dtype=complex) +start = time.clock() +for k in range(q,m): + AbatchDMD[:,:,k] = y[:,:k+1].dot(np.linalg.pinv(x[:,:k+1])) + evalsbatchDMD[:,k] = np.log(np.linalg.eigvals(AbatchDMD[:,:,k]))/dt +end = time.clock() +print "Batch DMD, time = " + str(end-start) + " secs" + + +# Online DMD, weighting = 1 +evalsonlineDMD1 = np.empty((n,m),dtype=complex) +odmd = dmdtools.OnlineDMD(n,1.0) +odmd.initialize(x[:,:q],y[:,:q]) +start = time.clock() +for k in range(q,m): + odmd.update(x[:,k],y[:,k]) + evalsonlineDMD1[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt +end = time.clock() +print "Online DMD, weighting = 1, time = " + str(end-start) + " secs" + + +# Online DMD, weighting = 0.9 +evalsonlineDMD2 = np.empty((n,m),dtype=complex) +odmd = dmdtools.OnlineDMD(n,0.9) +odmd.initialize(x[:,:q],y[:,:q]) +start = time.clock() +for k in range(q,m): + odmd.update(x[:,k],y[:,k]) + evalsonlineDMD2[:,k] = np.log(np.linalg.eigvals(odmd.A))/dt +end = time.clock() +print "Online DMD, weighting = 0.9, time = " + str(end-start) + " secs" + + +# visualize true, batch, online (weighting=1,0.9) +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsbatchDMD[0,q:]), 'r-',label='Batch',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD1[0,q:]), 'g--',label='Online, $wf=1$',linewidth=2.0) +plt.plot(t[q:], np.imag(evalsonlineDMD2[0,q:]), 'b--',label='Online, $wf=0.9$',linewidth=2.0) +plt.tick_params(labelsize=20) +plt.xlabel('Time', fontsize=20) +plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) +plt.legend(loc='best', fontsize=20, shadow=True) +plt.xlim([0,10]) +plt.ylim([1,2]) +plt.show() \ No newline at end of file diff --git a/python/scripts/window_demo.py b/python/scripts/window_demo.py new file mode 100644 index 0000000..3ea0b20 --- /dev/null +++ b/python/scripts/window_demo.py @@ -0,0 +1,135 @@ +# -*- coding: utf-8 -*- +""" +An example to demonstrate window dynamic mode decomposition + +We take a 2D time varying system given by dx/dt = A(t)x +where x = [x1,x2]', A(t) = [0,w(t);-w(t),0], +w(t)=1+epsilon*t, epsilon=0.1. The slowly time varying eigenvlaues of A(t) +are pure imaginary, +(1+0.1t)j and -(1+0.1t)j, where j is the imaginary unit. + +At time step k, define two matrix X(k) = [x(k-w+1),x(k-w+2),...,x(k)], +Y(k) = [y(k-w+1),y(k-w+2),...,y(k)], that contain the recent w snapshot pairs +from a finite time window, we would like to compute Ak = Yk*pinv(Xk). This can +be done by brute-force mini-batch DMD, and by efficient rank-2 updating window +DMD algrithm. For window DMD, at time k+1, we need to forget the old snapshot +pair xold = x(k-w+1), yold = y(k-w+1), and remember the new snapshot pair xnew += x(k+1), ynew = y(k+1). Mini-batch DMD computes DMD matrix by taking the +pseudo-inverse directly. Window DMD computes the DMD matrix by using efficient +rank-2 update idea. + +We compare the performance of window DMD with the brute-force mini-batch DMD +approach in terms of tracking time varying eigenvalues, by comparison with +the analytical solution. They should agree with each other (up to machine +round-offer errors). + +Authors: + Hao Zhang + Clarence W. Rowley + +References: + Hao Zhang, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta, + "Online Dynamic Mode Decomposition for Time-varying Systems," + arXiv preprint arXiv:1707.02876, 2017. + +Date created: April 2017 +""" + + +import sys +sys.path.append('..') + +import dmdtools +import numpy as np +from scipy.integrate import odeint +import time +import matplotlib.pyplot as plt + + +# define dynamics +epsilon = 1e-1 +def dyn(x,t): + x1, x2 = x + dxdt = [(1+epsilon*t)*x2,-(1+epsilon*t)*x1] + return dxdt +# integrate from initial condition [1,0] +tspan = np.linspace(0,10,101) +dt = 0.1 +x0 = [1,0] +xsol = odeint(dyn,x0,tspan).T +# extract snapshots +x, y = xsol[:,:-1], xsol[:,1:] +t = tspan[1:] +# true dynamics, true eigenvalues +n, m = len(x[:,0]), len(x[0,:]) +A = np.empty((n,n,m)) +evals = np.empty((n,m),dtype=complex) +for k in range(m): + A[:,:,k] = np.array([[0,(1+epsilon*t[k])],[-(1+epsilon*t[k]),0]]) + evals[:,k] = np.linalg.eigvals(A[:,:,k]) + + +# visualize snapshots +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(tspan, xsol[0,:], 'bs-', linewidth=2.0, label='$x_1(t)$') +plt.plot(tspan, xsol[1,:], 'g^-', linewidth=2.0, label='$x_2(t)$') +plt.legend(loc='best',fontsize=20 ,shadow=True) +plt.xlabel('Time', fontsize=20) +plt.title('Snapshots', fontsize=20) +plt.tick_params(labelsize=20) +plt.grid() +plt.show() + + +# mini-batch DMD, w = 10 +w = 10 +AminibatchDMD = np.empty((n,n,m)) +evalsminibatchDMD = np.empty((n,m),dtype=complex) +start = time.clock() +for k in range(w,m): + AminibatchDMD[:,:,k] = y[:,k-w+1:k+1].dot(np.linalg.pinv(x[:,k-w+1:k+1])) + evalsminibatchDMD[:,k] = np.log(np.linalg.eigvals(AminibatchDMD[:,:,k]))/dt +end = time.clock() +print "Mini-batch DMD, w = 10, time = " + str(end-start) + " secs" + + +# Window DMD, w = 10, weighting = 1 +evalswindowDMD1 = np.empty((n,m),dtype=complex) +wdmd = dmdtools.WindowDMD(n,w,1) +wdmd.initialize(x[:,:w],y[:,:w]) +start = time.clock() +for k in range(w,m): + wdmd.update(x[:,k],y[:,k]) + evalswindowDMD1[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt +end = time.clock() +print "Window DMD, w=10, weighting = 1, time = " + str(end-start) + " secs" + + +# Window DMD, w = 10, weighting = 0.5 +evalswindowDMD2 = np.empty((n,m),dtype=complex) +wdmd = dmdtools.WindowDMD(n,w,0.5) +wdmd.initialize(x[:,:w],y[:,:w]) +start = time.clock() +for k in range(w,m): + wdmd.update(x[:,k],y[:,k]) + evalswindowDMD2[:,k] = np.log(np.linalg.eigvals(wdmd.A))/dt +end = time.clock() +print "Window DMD, w=10, weighting = 0.5, time = " + str(end-start) + " secs" + + +# visualize true, batch, window +plt.figure() +plt.rc('text', usetex=True) +plt.rc('font', family='serif') +plt.plot(t, np.imag(evals[0,:]), 'k-',label='True',linewidth=2.0) +plt.plot(t[w:], np.imag(evalsminibatchDMD[0,w:]), 'r-',label='Mini-batch, $w=10$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD1[0,w:]), 'g--',label='Window, $w=10$, $wf=1$',linewidth=2.0) +plt.plot(t[w:], np.imag(evalswindowDMD2[0,w:]), 'b--',label='Window, $w=10$, $wf=0.5$',linewidth=2.0) +plt.tick_params(labelsize=20) +plt.xlabel('Time', fontsize=20) +plt.ylabel('Im($\lambda_{DMD}$)', fontsize=20) +plt.legend(loc='best', fontsize=20, shadow=True) +plt.xlim([0,10]) +plt.ylim([1,2]) +plt.show() \ No newline at end of file