emgr v5.1 release

CITATION

@@ -1,13 +1,13 @@
 Cite As:
 
-C. Himpe (2016). emgr - Empirical Gramian Framework (Version 5.0) [Software]. Available from http://gramian.de . doi: 10.5281/zenodo.162135 .
+C. Himpe (2017). emgr - EMpirical GRamian Framework (Version 5.1) [Software]. Available from http://gramian.de . doi:10.5281/zenodo.580804
 
 BibTeX Entry:
 
 @MISC{emgr,
  author = {C. Himpe},
- title = {{emgr - Empirical Gramian Framework (Version 5.0)}},
+ title = {{emgr - EMpirical GRamian Framework (Version 5.1)}},
  howpublished = {\url{http://gramian.de}},
- year = {2016},
- doi = {10.5281/zenodo.162135}
+ year = {2017},
+ doi = {10.5281/zenodo.580804}
 } 

CODE

@@ -1,8 +1,8 @@
 name: Empirical Gramian Framework
 shortname: emgr
-version: 5.0
-release-date: 2016-10-20
-id: 10.5281/zenodo.162135
+version: 5.1
+release-date: 2017-05-18
+id: 10.5281/zenodo.580804
 id-type: doi
 author: Christian Himpe
 orcid: 0000-0003-2194-6754
@@ -13,7 +13,8 @@ license-type: open
 repository: github.com/gramian/emgr
 repository-type: git
 language: Matlab
-dependencies: Octave >=4.0, Matlab >=2016b
+dependencies: Octave >=4.2, Matlab >=2016b
 systems: Linux, Windows
 website: gramian.de
 keywords: Controllability, Observability, Reduced Order Modelling, Model Reduction, Model Order Reduction
+msc: 93b99

LICENSE

@@ -1,6 +1,6 @@
 All source code is licensed under the open source BSD 2-clause license:
 
-Copyright (c) 2013-2016, Christian Himpe
+Copyright (c) 2013--2017, Christian Himpe
 
 All rights reserved.
 

README.md

@@ -1,15 +1,132 @@
-![emgr Logo](emgr.png) emgr
-===========================
-Empirical Gramian Framework (Version 5.0) ![DOI badge](https://zenodo.org/badge/doi/10.5281/zenodo.162135.png)
+![emgr logo](emgr.png) emgr - EMpirical GRamian Framework
+=========================================================
 
-Empirical gramians can be computed for linear and nonlinear state-space control systems for purposes of model order reduction (MOR), system identification (SYSID) and uncertainty quantification (UQ).
-Model reduction using empirical gramians can be applied to the state-space, to the parameter-space or both through combined reduction.
-For state reduction, balanced truncation of the empirical controllability gramian and the empirical observability gramian, or alternatively, direct truncation (approximate balancing) of the empirical cross gramian (or the empirical linear cross gramian for large-scale linear systems) is available.
-For parameter reduction, parameter identification and sensitivity analysis the empirical sensitivity gramian (controllability of parameters) or the empirical identifiability gramian (observability of parameters) are provided.
-Combined state and parameter reduction is enabled by the empirical joint gramian, which computes controllability and observability of states (cross gramian) and observability of parameters (cross-identifiability gramian) concurrently.
-The empirical gramian framework - emgr is a compact open-source toolbox for (empirical) GRAMIAN-based model reduction and compatible with OCTAVE and MATLAB.
-emgr provides a common interface for the computation of empirical gramians and empirical covariance matrices.
+![code meta-data.](https://img.shields.io/badge/code_meta--data-%E2%9C%93-brightgreen.svg) 
+![zenodo listed.](https://zenodo.org/badge/doi/10.5281/zenodo.580804.png)
+![matlab compatible](https://img.shields.io/badge/matlab-compatible-lightgrey.svg)
 
-More information at: http://gramian.de
+* emgr - EMpirical GRamian Framework ( gramian.de )
+* version: **5.1** ( 2017-05-18 )
+* by: Christian Himpe ( 0000-0003-2194-6754 )
+* under: BSD 2-Clause License ( opensource.org/licenses/BSD-2-Clause )
+* summary: Empirical Gramians for model reduction of input-output systems.
 
-Code Meta Information: [CODE](CODE)
+## Scope
+
+* Model Reduction / Model Order Reduction (MOR)
+  * parametric Model Order Reduction (pMOR) / Robust Reduction
+  * nonlinear Model Order Reduction (nMOR)
+  * Parameter Identification / Parameter Reduction
+  * **Combined State and Parameter Reduction** (Combined Reduction)
+* Sensitivity Analysis
+* Decentralized Control
+* Optimal Sensor Placement / Optimal Actuator Placement
+
+## Empirical Gramians (Empirical Gramian Matrix Types)
+
+* Empirical Controllability Gramian
+* Empirical Observability Gramian
+* Empirical Cross Gramian
+* Empirical Linear Cross Gramian
+* Empirical Sensitivity Gramian
+* Empirical Identifiability Gramian
+* Empirical Joint Gramian (Empirical Cross-Identifiability Gramian)
+
+### Features
+
+* Interfaces for:
+  * custom solvers / integrators
+  * custom inner products / dot products
+  * distributed memory / column-wise computation
+* Non-symmetric cross gramian option for:
+  * Empirical Cross Gramian
+  * Empirical Linear Cross Gramian
+  * Empirical Joint Gramian
+* Compatible with:
+  * GNU Octave
+  * Mathworks MATLAB
+* Vectorized and parallelizable
+
+## Basic Usage
+
+Run a minimal example in a Matlab interpreter like OCTAVE or MATLAB:
+```
+RUNME
+```
+
+To use specific utilities, tests or demos include:
+```
+addpath('utils')
+addpath('tests')
+addpath('demos')
+```
+
+## Files and Folders
+
+[`README.md`](README.md) Basic Information
+
+[`RUNME.m`](RUNME.m) Minimal Code Example
+
+[`CODE`](CODE) Code Meta Information
+
+[`CITATION`](CITATION) Citation Information
+
+[`LICENSE`](LICENSE) License Indormation
+
+[`emgr.m`](emgr.m) Empirical Gramian Framework (Main File)
+
+[`emgr_oct.m`](emgr_oct.m) Empirical Gramian Framework (Optional Octave Variant) 
+
+[`emgr_lgc.m`](emgr_lgc.m) Empirical Gramian Framework (Pre 2016b Matlab)
+
+[`emgrtest.m`](emgrtest.m) Sanity test via experimental order of convergence
+
+[`testall.m`](testall.m) Run all tests and demos
+
+[`emgr_ref.pdf`](emgr_ref.pdf) emgr cheat sheet
+
+[`utils`](utils) Folder containing utilities
+
+[`tests`](tests) Folder containing tests
+
+[`demos`](demos) Folder containing demos
+
+## Documentation ToC
+
+### [Summary](http://gramian.de/#summary)
+
+### [Scope](http://gramian.de/#scope)
+
+### [Download](http://gramian.de/#download)
+
+### [License](http://gramian.de/#license)
+
+### [Disclaimer](http://gramian.de/#disclaimer)
+
+### [Usage](http://gramian.de/#usage)
+
+### [Arguments](http://gramian.de/#arguments)
+
+### [Empirical Gramian Types](http://gramian.de/#gramians)
+
+### [Option Flags](http://gramian.de/#options)
+
+### [Custom Solver](http://gramian.de/#solver)
+
+### [Extra Utilities](http://gramian.de/#utilities)
+
+### [Tests](http://gramian.de/#tests)
+
+### [Demos](http://gramian.de/#demos)
+
+### [About](http://gramian.de/#about)
+
+### [References](http://gramian.de/#references)
+
+### [Links](http://gramian.de/#links)
+
+### [Notes](http://gramian.de/#notes)
+
+### [Troubleshooting](http://gramian.de/#troubleshooting)
+
+## More information at: http://gramian.de

RUNME.m

@@ -0,0 +1,22 @@
+%%% summary: RUNME (Minimal setup and test script)
+%%% project: emgr - Empirical Gramian Framework ( gramian.de )
+%%% authors: Christian Himpe ( 0000-0003-2194-6754 )
+%%% license: 2-Clause BSD (2017)
+%$
+
+% Linear System
+A = -eye(4)
+B = [0;1;0;1]
+C = [0,0,1,1]
+
+% (Empircal) Controllability Gramian
+emgr(@(x,u,p,t) A*x+B*u,@(x,u,p,t) C*x,[1,4,1],[0.1,1.0],'c')
+
+% (Empirical) Observability Gramian
+emgr(@(x,u,p,t) A*x+B*u,@(x,u,p,t) C*x,[1,4,1],[0.1,1.0],'o')
+
+% (Empirical) Cross Gramian
+emgr(@(x,u,p,t) A*x+B*u,@(x,u,p,t) C*x,[1,4,1],[0.1,1.0],'x')
+
+% (Empirical) Linear Cross Gramian
+emgr(@(x,u,p,t) A*x+B*u,@(x,u,p,t) A'*x+C'*u,[1,4,1],[0.1,1.0],'y')

demos/advection.m

@@ -2,7 +2,7 @@ function advection(o)
 %%% summary: advection (finite difference discretized transport equation)
 %%% project: emgr - Empirical Gramian Framework ( gramian.de )
 %%% authors: Christian Himpe ( 0000-0003-2194-6754 )
-%%% license: 2-Clause BSD (2013--2016)
+%%% license: 2-Clause BSD (2013--2017)
 %$
     if(exist('emgr')~=2)
         error('emgr not found! Get emgr at: http://gramian.de');
@@ -12,42 +12,43 @@ function advection(o)
         fprintf('emgr (version: %1.1f)\n',emgr('version'));
     end
 
-%% SETUP
-    M = 0;
-    N = 256;
-    Q = N;
-    R = 10;
-    T = [0.01,0.1];
-    L = floor(T(2)/T(1)) + 1;
-    U = @(t) 0;
-    X = exp(-linspace(-2,8,N).^2)';
+%% SYSTEM SETUP
+    M = 0;				% number of inputs
+    N = 256;				% number of states
+    Q = N;				% number of outputs
+    R = 10;				% target reduced order
+    h = 0.01;				% time step size
+    T = 0.1;				% time horizon
+    X = exp(-linspace(-2,8,N).^2)';	% initial state
+    U = @(t) 0;				% zero input function
+    P = 0.55;				% parameter
 
-    P = 0.55;
-    A = spdiags(N*[ones(N,1),-ones(N,1)],[-1,0],N,N);
+    A = spdiags(N*[ones(N,1),-ones(N,1)],[-1,0],N,N);	% system matrix
 
-    LIN = @(x,u,p,t) p*A*x;
-    ADJ = @(x,u,p,t) p*A'*x + u;
-    OUT = @(x,u,p,t) x;
+    LIN = @(x,u,p,t) p*A*x;		% vector field
+    ADJ = @(x,u,p,t) p*A'*x + u;	% adjoint vector field
+    OUT = @(x,u,p,t) x;			% output functional
 
-%% OFFLINE
+%% REDUCED ORDER MODEL PROJECTION ASSEMBLY
     tic;
-    WO = emgr(ADJ,1,[N,N,Q],T,'c',P);
-    [UU,D,VV] = svd(WO); UU = UU(:,1:R); VV = UU';
+    WO = emgr(ADJ,1,[N,N,Q],[h,T],'c',P);
+    [UU,D,VV] = svd(WO);
+    UU = UU(:,1:R);
+    VV = UU';
     a = VV*A*UU;
     c = UU;
     x = VV*X;
     lin = @(x,u,p,t) p*a*x;
     out = @(x,u,p,t) c*x;
     OFFLINE = toc
 
-%% ONLINE
+%% REDUCED ORDER MODEL SIMULATION
     y = ODE(lin,out,[0.01,1.0],x,U,P);
 
-%% OUTPUT
+%% PLOT REDUCED ORDER MODEL PHASE SPACE
     if(nargin>0 && o==0), return; end; 
     figure('Name',mfilename,'NumberTitle','off');
     imagesc(sparse(y)); caxis([0,max(y(:))]);
-    if(exist('viridis')==0), colormap(hot); end;
     set(gca,'YTick',0,'xtick',[]); ylabel('X'); xlabel('t');
     pbaspect([2,1,1]);
     if(nargin>0 && o==1), print('-dpng',[mfilename(),'.png']); end;

demos/benchmark_fss.m

@@ -0,0 +1,91 @@
+function benchmark_fss(o)
+%%% summary: benchmark_fss (Flexible Space Structures Benchmark)
+%%% project: emgr - Empirical Gramian Framework ( http://gramian.de )
+%%% authors: Christian Himpe ( 0000-0003-2194-6754 )
+%%% license: 2-Clause BSD (2017)
+%$
+    if(exist('emgr')~=2)
+        error('emgr not found! Get emgr at: http://gramian.de');
+    else
+        global ODE;
+        ODE = [];
+        fprintf('emgr (version: %1.1f)\n',emgr('version'));
+    end
+
+    rand('seed',1009);
+
+%% SETUPs
+
+    K = 32;				% number of modes
+    xi = rand(1,K)*0.001;		% damping ratio
+    omega = rand(1,K)*100;		% natural frequencies
+
+    N = 2*K;				% number of states
+    M = 1;				% number of inputs
+    Q = M;				% number of outputs
+
+    A_k = cellfun(@(p) sparse([-2.0*p(1)*p(2),-p(2);p(2),0]), ...
+                  num2cell([xi;omega],1),'UniformOutput',0);
+    A = blkdiag(A_k{:});		% system matrix
+    B = kron(rand(K,M),[1;0]);		% input matrix
+    C = 10.0*rand(Q,2*K);		% output matrix
+
+    h = 0.001;				% time step size
+    T = 1.0;				% time horizon			
+    L = floor(T/h) + 1;			% number of time steps
+    U = @(t) ones(M,1)*(t<=h)/h;	% impulse input function
+    X = zeros(N,1);
+
+    LIN = @(x,u,p,t) A*x + B*u;		% vector field
+    OUT = @(x,u,p,t) C*x;		% output functional
+
+%% FULL ORDER
+    Y = ODE(LIN,OUT,[h,T],X,U,0);
+    %figure; plot(0:h:T,Y); return;
+    n1 = norm(Y(:),1);
+    n2 = norm(Y(:),2);
+    n8 = norm(Y(:),Inf);
+
+%% OFFLINE
+    tic;
+    s = 16;
+    P = ceil(N/s);
+    wx = cell(1,Q);
+    for p=1:P
+        wx{p} = emgr(LIN,OUT,[M,N,Q],[h,T],'x',0,[0,0,0,0,0,0,0,0,0,0,s,p]);
+    end;
+    [UU,DD,VV] = svd(cell2mat(wx));
+    OFFLINE_TIME = toc
+
+%% EVALUATION
+    l1 = zeros(1,N-1);
+    l2 = zeros(1,N-1);
+    l8 = zeros(1,N-1);
+
+    for n=1:N-1
+        uu = UU(:,1:n);
+        a = uu'*A*uu;
+        b = uu'*B;
+        c = C*uu;
+        x = uu'*X;
+        lin = @(x,u,p,t) a*x + b*u;
+        out = @(x,u,p,t) c*x;
+        y = ODE(lin,out,[h,T],x,U,0);
+        l1(n) = norm(Y(:)-y(:),1)/n1;
+        l2(n) = norm(Y(:)-y(:),2)/n2;
+        l8(n) = norm(Y(:)-y(:),Inf)/n8;
+    end;
+
+%% OUTPUT
+    if(nargin>0 && o==0), return; end; 
+    figure('Name',mfilename,'NumberTitle','off');
+    semilogy(1:N-1,l1,'r','linewidth',2); hold on;
+    semilogy(1:N-1,l2,'g','linewidth',2);
+    semilogy(1:N-1,l8,'b','linewidth',2); hold off;
+    xlim([1,N-1]);
+    ylim([1e-16,1]);
+    pbaspect([2,1,1]);
+    legend('L1 Error ','L2 Error ','L8 Error ','location','northeast');
+    if(nargin>0 && o==1), print('-dsvg',[mfilename(),'.svg']); end;
+end
+