Keedmd continuous

scripts/Koopman-ID/DMD/CompanionMatrix_DMD.m

@@ -0,0 +1,54 @@
+function [ KModes,KEv,Norms ] = CompanionMatrix_DMD( Data )
+% Dynamic Mode Decomposition as presented by 
+% "Spectral analysis of nonlinear flows" by Rowley et al., Journal of FLuid
+% Mechanics, 2009
+
+
+% inputs : 
+% Data - the data matrix: each column is a a set of measurements done at
+% each instant - the sampling frequency is assumed to be constant
+
+% outputs: 
+% 1 - KModes- Koopman or Dynamic Modes: structures that scale exponentially
+% with time - the exponent is given by Koopman or Dynamic Eigenvalues and
+% could be imaginary as well
+
+% 2- KEv - Koopman or Dynamic Eigenvalues: the exponents and frequencies
+% that make up the time-varaiation of data in time
+
+% 3- Norms - Euclidean (vector) norm of each mode, used to sort the data
+
+
+
+
+X=Data(:,1:end-1);
+
+c = pinv(X)*Data(:,end);
+
+m = size(Data,2)-1;
+C = spdiags(ones(m,1),-1,m,m);
+C(:,end)=c;                  %% companion matrix
+
+[P,D]=eig(full(C));                           %% Ritz evalue and evectors
+
+KEunsrtd = diag(D);                     %% Koopman Evalues
+KMunsrtd = X*P;                         %% Koopman Modes
+
+Cnorm = sqrt(sum(abs(KMunsrtd).^2,1));  %% Euclidean norm of Koopman Modes
+[Norms, ind]=sort(Cnorm,'descend');      %% sorting based on the norm
+
+KEv= KEunsrtd(ind);           %% Koopman Eigenvalues
+
+KModes = KMunsrtd(:,ind);                %% Koopman Modes
+
+
+
+end
+
+
+%=========================================================================%
+% Hassan Arbabi - 08-17-2015
+% Mezic research group
+% UC Santa Barbara
+% arbabiha@gmail.com
+%=========================================================================%

scripts/Koopman-ID/DMD/Exact_DMD.m

@@ -0,0 +1,71 @@
+function [ ProjectedModes,DEv,ExactModes,Norm ] = Exact_DMD( X,Y,varargin )
+% Dynamic Mode Decomposition as presented by 
+% "On Dynamic Mode Decomposition: theory and applications" by Tu et al.,
+% arXiv, 2013
+
+
+
+% inputs : 
+% Data Sets X and Y- should have the same size 
+% each column of X is a set of measurements done at an instant 
+% each column of Y is the image of the corresponding columns from X
+
+% or 
+
+% ( X,Y,Tol ) - with Tol (optional) being the threshold for filtering thru SVD - the
+% default value is 1e-10
+
+
+
+% outputs: 
+% 1 - Projected Dynamic Modes
+% 2 - Exact Dynamic Modes
+% 3 - Dynamic Eigenvalues
+% 4 - Norms - Euclidean (vector) norm of each mode, used to sort the data
+
+% setting SVD hard threshold
+if isempty(varargin)
+    Tol=1e-10;
+else
+    Tol=varargin{1};
+end
+
+
+disp(['Tolerance used for filtering in DMD:',num2str(Tol)])
+
+[U,S,V]=svd(X,'econ');
+
+k = find(diag(S)>Tol,1,'last');
+disp(['DMD subspace dimension:',num2str(k)])
+
+U = U(:,1:k); V=V(:,1:k); S=S(1:k,1:k);
+
+Atilde = ((U'*Y) *V )* diag((1./diag(S)));
+
+
+[W,DEv]=eig(Atilde);
+DEv = diag(DEv);
+
+ProjectedModes = U*W;
+
+ExactModes = bsxfun(@times,1./DEv.',((Y*V)*S^(-1))*W);
+
+
+
+if nargout>3
+   b = pinv(ProjectedModes)*X(:,end);     % terminal coordinates in the Koopman subspace
+    [Norm,Index]=sort(abs(b),'descend');
+    DEv = DEv(Index);
+    ProjectedModes = ProjectedModes(:,Index);
+    disp('modes sorted based on energy contribution to the last snapshot')
+end
+NORM=Norm
+end
+
+
+%=========================================================================%
+% Hassan Arbabi - 08-17-2015
+% Mezic research group
+% UC Santa Barbara
+% arbabiha@gmail.com
+%=========================================================================%

scripts/Koopman-ID/DMD/Hankel_DMD.m

@@ -0,0 +1,65 @@
+function [ HModes, HEvalues, Norms ] = Hankel_DMD( Data,n,m,varargin )
+%HANKEL-Dynamic Mode Decompsition as presented in
+% "Ergodic theory, dynamic mode decomposition and computation of Koopman
+% spectral properties" by H. Arbabi and I. Mezic, arXiv:1611.06664
+
+% inputs : 
+% Data - the data matrix: each row is a time series data on an observable
+% m,n: number of columns and rows, respectively, in the data Hankel
+% matrices: size(Data,2)>m+n and preferably m>>n 
+
+% or 
+
+% ( Data,n,m,tol )- with Tol (optional) being the threshold for filtering thru SVD - the
+% default value is 1e-10
+
+
+
+
+% outputs: 
+% 1 - HModes- Dynamic modes which approximate the Koopman eigenfunctions
+
+% 2- HEvalues - Koopman or Dynamic Eigenvalues: the exponents and frequencies
+% that make up the time-varaiation of data in time
+
+% 3- Norms - The L2-norm of Koopman eigenfunction contribution to the
+% observable
+
+
+% setting the SVD hard threshold for DMD
+if isempty(varargin)
+    Tol=1e-10;
+else
+    Tol=varargin{1};
+end
+
+index1 = 1:n;
+index2 = n:n+m-1;
+
+X = []; Y=[];
+
+for ir = 1:size(Data,1)
+
+    % Hankel blocks ()
+    c = Data(ir,index1).'; r = Data(ir,index2);
+    H = hankel(c,r).';
+    c = Data(ir,index1+1).'; r = Data(ir,index2+1);
+    UH= hankel(c,r).';
+
+    % scaling of Hankel blocks
+    if ir>1
+        alpha = norm(H(:,1))/norm(X(:,1));
+        H = alpha*H;
+        UH= alpha* UH;
+    end
+
+    % the data matrices fed to exact DMD
+    X=[X,H]; Y=[Y,UH];
+
+end
+    % the DMD
+   [HModes,HEvalues,~,Norms ] = DMD.Exact_DMD( X,Y,Tol );
+
+
+end
+

scripts/Koopman-ID/DMD/HarmonicAverage.m

@@ -0,0 +1,22 @@
+function [ Average ] = HarmonicAverage( Data,w,t,Filter)
+% computes the harmonic average using the Filter
+
+HarmonicWeight = 2*exp(-1i*w'*t);            
+
+% using the specified filter
+if exist('Filter','var')
+   switch Filter
+       case 'Hamming'
+           HarmonicWeight=bsxfun(@times,HarmonicWeight,hamming(length(t))');
+       case 'Hann'
+           HarmonicWeight=bsxfun(@times,HarmonicWeight,hann(length(t))');
+       case 'exponential'       % the exponential weighting defined by Das & Yorke 2015
+           HarmonicWeight=bsxfun(@times,HarmonicWeight,ExpWeight(length(t))*length(t));
+       otherwise
+           disp('filter not defined ... no filter used')
+   end 
+end
+
+Average = HarmonicWeight*Data'./length(t);     
+
+end

scripts/Koopman-ID/DMD/SVDenhanced_DMD.m

@@ -0,0 +1,74 @@
+function [ DModes,DEv,Norm ] = SVDenhanced_DMD( Data,varargin )
+% Dynamic Mode Decomposition as presented by 
+% "Dynamic Mode Decomposition of numerical and experimental data" by Peter
+% J. Schmid, Journal of Fluid Mechanics, 2010
+
+% inputs : 
+% Data - the data matrix: each column is a a set of measurements done at
+% each instant - the sampling frequency is assumed to be constant
+
+% or 
+
+% ( Data,Tol ) - with Tol (optional) being the threshold for filtering thru SVD - the
+% default value is 1e-10
+
+% outputs: 
+% 1 - DModes- Koopman or Dynamic Modes: structures that scale exponentially
+% with time - the exponent is given by Koopman or Dynamic Eigenvalues and
+% could be imaginary as well
+
+% 2- DEv - Koopman or Dynamic Eigenvalues: the exponents and frequencies
+% that make up the time-varaiation of data in time
+
+% 3- Norm - The norm here is set to the energy contribution of each mode to
+% the last snapshot of data and then used to sort the modes
+
+
+
+
+tic
+if isempty(varargin)
+    Tol=1e-10;
+else
+    Tol=varargin{1};
+end
+
+disp(['Tolerance used for filtering in DMD:',num2str(Tol)])
+
+X=Data(:,1:end-1);
+Y=Data(:,2:end);
+
+
+[U,S,V]=svd(X,0);
+r = find(diag(S)>Tol,1,'last');
+disp(['DMD subspace dimension:',num2str(r)])
+
+U = U(:,1:r);
+S = S(1:r,1:r);
+V = V(:,1:r);
+
+Atilde = (U'*(Y*V)) / S; 
+
+[w,lambda]=eig(Atilde);
+
+DModes = U*w;
+DEv = diag(lambda);
+
+if nargout>2
+    b = pinv(DModes)*Data(:,end);     % terminal coordinates in the Koopman subspace
+    [Norm,Index]=sort(abs(b),'descend');
+    DEv = DEv(Index);
+    DModes = DModes(:,Index);
+    disp('modes sorted based on energy contribution to the last snapshot')
+end
+
+disp(['time to compute DMD:',num2str(toc)])
+
+end
+
+%=========================================================================%
+% Hassan Arbabi - 08-17-2015
+% Mezic research group
+% UC Santa Barbara
+% arbabiha@gmail.com
+%=========================================================================%

scripts/Koopman-ID/KEEDMD.m

@@ -0,0 +1,82 @@
+function [A_koop, B_koop, C_koop, liftFun] = KEEDMD(n,m,Ntraj, Ntime, N_basis,...
+        pow_eig_pairs, Xstr, Ustr, qf, A_nom, B_nom, K_nom, deltaT)
+    %Identify lifted state space model using Extended Dynamic Mode
+    %Decomposition
+
+    %Inputs:
+    %   n               - Number of states
+    %   m               - Number of control inputs
+    %   N_basis         - Number of basis functions to use
+    %   basis_function  - Type of basis functions (only rbfs implemented)
+    %   rbf-type        - Type of rbf function 
+    %   center_type     - What centers to use for rbfs ('data'/'random') 
+    %   eps             - Width of rbf if rbf type has width parameter
+    %   plot_basis      - (true/false) determines to plot basis functions
+    %   xlim            - Plot limits x-axis
+    %   ylim            - Plot limits y-axis
+    %   Xacc            - Trajectory data
+    %   Yacc            - Trajectory data at next time step
+    %   Uacc            - Control input for trajectory data
+
+    %Outputs:
+    %   A_edmd          - Passive dynamics matrix in lifted space
+    %   B_edmd          - Actuation matrix in lifted space
+    %   C_edmd          - Projection matrix from lifted space to outputs
+
+    % ************************** Prepare data *****************************
+    Xstr_shift = zeros(size(Xstr)); %Shift dynamics such that origin is fixed point
+    X = [];
+    X_dot = [];
+    U = [];
+    for i = 1 : Ntraj
+       Xstr_shift(:,i,:) = Xstr(:,i,:)-qf(:,i);
+       X = [X reshape(Xstr_shift(:,i,:),size(Xstr_shift,1),size(Xstr_shift,3))];
+       X_dot = [X_dot num_diff(reshape(Xstr_shift(:,i,:),size(Xstr_shift,1),size(Xstr_shift,3)),deltaT)];
+       U = [U reshape(Ustr(:,i,:),size(Ustr,1),size(Ustr,3))];
+    end
+
+    % ******************** Construct eigenfunctions ***********************
+    [A_eigfuncs, liftFun] = construct_eigfuncs(n, N_basis,pow_eig_pairs, ...
+                                            A_nom, B_nom, K_nom, X, X_dot);
+
+    %Perform lifting
+    liftFun = @(xx) [xx; liftFun(xx)];
+    Nlift = length(liftFun(zeros(n,1)));
+    Xlift = [];
+    Xlift_dot = [];
+    for i = 1 : Ntraj
+       Xlift_temp = liftFun(reshape(Xstr_shift(:,i,:),size(Xstr_shift,1),size(Xstr_shift,3)));
+       Xlift = [Xlift Xlift_temp];
+       Xlift_dot = [Xlift_dot num_diff(Xlift_temp,deltaT)];
+    end
+
+    % ********************** Build predictor *********************************
+
+    %Set up regression for A and B:
+    X_vel = [Xlift; U];
+    Y_vel = Xlift_dot(n/2+1:n,:);
+    A_vel_x = lasso(X_vel',Y_vel(1,:)','Lambda',1e-3, 'Alpha', 0.5);
+    A_vel_th = lasso(X_vel',Y_vel(2,:)','Lambda',1e-3, 'Alpha', 0.5);
+
+    %Perform regression and enforce known structure:
+    A_koop = zeros(Nlift);
+    A_koop(1:n/2,:) = [zeros(n/2) eye(n/2) zeros(n/2,size(A_koop,2)-n)];
+    A_koop(n/2+1:n,:) = [A_vel_x(1:Nlift,:)'; A_vel_th(1:Nlift,:)']; 
+    A_koop(n+1:end,n+1:end) = A_eigfuncs;
+
+    Y_kin = X_dot(1:n/2,:) - A_koop(1:n/2,:)*Xlift;
+    X_kin = U;
+    B_kin = X_kin'\Y_kin';
+
+    Y_eig = Xlift_dot(n+1:end,:) - A_eigfuncs*Xlift(n+1:end,:);
+    X_eig = U-K_nom*X;
+    B_eig = X_eig'\Y_eig';
+
+    B_koop = [B_kin'; A_vel_x(Nlift+1:end,:); A_vel_th(Nlift+1:end,:); B_eig'];  
+
+    C_koop = zeros(n,size(A_koop,1));
+    C_koop(1:n,1:n) = eye(n);
+
+    %Subtract effect of nominal controller from A:
+    A_koop(n+1:end,1:n) = A_koop(n+1:end,1:n)-B_eig'*K_nom;
+end