MOR基础笔记

Krylov Subspace

使用幂迭代查找A^{m\times m}的最大特征值\lambda_1的一种方法是：\forall b_0，\lim_{n\to+\infin}A^nb_0=Av_1，v_1是\lambda_1的特征向量。

推广到一个Krylov矩阵，

K_n=[b\quad Ab\quad A^2b\quad\cdots\quad A^{n-1}b],

我们期望正交后，得到的向量集是前n大特征值对应特征向量的良好近似。

Arnordi Iteration

Arnordi迭代使用改进的GS正交化方法达成这正交化Krylov子空间的目的（这个h下标写得有一定误导性其实，比如h_{1,1}实际上是向量1和2的内积）：

\begin{aligned} &\textbf{Start from}\ q_1\ \text{with norm 1, for k from 2 to n.}\\ &\quad q_k=Aq_{k-1}\\ &\text{Then do Gram-Schmidt:}\\ &\quad\textbf{for j from 1 to k − 1:}\\ &\qquad h_{j,k-1}:=\langle q_j,q_k\rangle\\ &\qquad q_k=q_k-h_{j,k-1}q_j\\ &\quad h_{k,k-1}:=\langle q_k,q_k\rangle\\ &\quad q_k=q_k/h_{k,k-1} \end{aligned}

将h按照下标排列起来你会发现能得到一个上海森堡矩阵H_k=\{h_{i,j}\}，其中（思考一下AQ乘完之后是什么就比较显然了）：

H_k=Q_k^*AQ_k,\quad Q_k:=[q_1\quad \cdots\quad q_n].\\

如果继续迭代：

AQ_k=Q_{k+1}\tilde{H_k},\quad \tilde{H_k}=[H_k^T\quad[0 \cdots h_{k+1,k}]^T]^T

拓展

如果A是一个Hermitian（复共轭对称）矩阵，那么Arnordi算法变成Lanczos算法。
GMRES是基于Arnordi算法的线性求解算法。

Block Arnordi

MOR

morwiki：https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Main_Page

默认的，本笔记研究的降阶形式为：

\begin{aligned} \boldsymbol{C}\dot{x_n}+\boldsymbol{G}x_n&=\boldsymbol{B}u_p\\ i_p&=\boldsymbol{L^T}x_n \end{aligned} \tag{MOR}

n是节点支路规模，p是端口规模，一般地，如果需要有效降阶，默认n>>p.

PRIMA

我现在对于PRIMA的创新点理解是这样子的，首先，利用Krylov子空间，可以得到

\begin{aligned} \boldsymbol{A}:=\boldsymbol{-G^{-1}C}&\quad\boldsymbol{R}:=\boldsymbol{G^{-1}B}\\ colsp(\boldsymbol{X})&=kr(\boldsymbol{A},\boldsymbol{R},q)\\ \boldsymbol{X^TAX}&=\boldsymbol{H_q}\\ \boldsymbol{X^TX}&=\boldsymbol{I_q} \end{aligned}

正常来说的约简方式是找一个x_n=\boldsymbol{X}z_q：

\begin{aligned} -\boldsymbol{X^TG^{-1}CX}\dot{z_q}&=\boldsymbol{X^TX}z_q-\boldsymbol{X^TG^{-1}B}u_p\\ i_p&=\boldsymbol{L^TX}z_q \end{aligned}

但是PRIMA不，它不用这种可推导的方式，他直接：

\begin{aligned} -\boldsymbol{X^TCX}\dot{z_q}&=\boldsymbol{X^TGX}z_q-\boldsymbol{X^TB}u_p\\ i_p&=\boldsymbol{L^TX}z_q \end{aligned}

然后它再去证明矩匹配和passivity（暂时没看详细）那一堆东西，我现在看PRIMA启发式约简也好，直接推导式约简也好，有一个问题就是q<n，他把这个x给变元了，我感觉DeMOR做的一件事情就是保留信息，它只去看一个端口，那他的n=1，q是多少并不是很重要了，反正都是用其他的信息去生成一个端口的内容，信息损失可能就没有那么大。

破案了，真正的原因是DeMOR由于实际上保留的端口少，那么在同样的\boldsymbol{X}的列数下，DeMOR保留的矩更多。

RGA

假设H(s)是传输函数，Y(s)=H(s)U(s)，那么RGA定义为：

\Lambda(H)=\{\lambda_{ij}\}=\{\frac{\partial y_i/\partial u_j}{\partial u_j/\partial y_i}\}=H(0)\circ H(0)^{-T}

这个求解需要求逆，代价非常大。

DeMOR

本质上是利用RGA做解耦，将一定阈值以上的端口作为主要端口提取出来（s for selected）：

B_i:=[b_{s1}\quad \cdots\quad b_{sm}]\qquad L_i:=l_{\ \text{observed}}

那么对于某个给定的端口i，传输函数就被简化成

H_i(s)=L_i^T(sC+G)^{-1}B_i\tag{DeMOR a}

只需要对这个解耦的传输方程做PRIMA降阶，就可以得到被观察端口的降阶形式，这个时候我们观察解耦子空间colsp(\boldsymbol{X})=kr(\boldsymbol{A},\boldsymbol{R_i},q)，会发现由于\boldsymbol{R_i}=G^{-1}B_i，子空间的维度从q\times p，变成了q\times m，m是B的列数，且在期望上至少是p的十分之一以下——郝博指出了这一点，这时候翻看代码就会发现论文中PRIMA的效果不如DeMOR的原因：在相同的子空间列数的前提下，DeMOR匹配到的矩就是更多。

但是没有任何数学证明这个问题之外，本文仍然有一个巨大问题：它用B_i 得到的\tilde B究竟在多大程度上拟合了B？这是不知道的。最离谱的地方在于，它最后降阶的那个传输函数并非式（DeMOR a），就是说，实际上降阶后它的\tilde B:=X_i^TB\neq X_i^TB_i，也就是子空间明明没有用到未被挑选出来的端口的信息，但最后却用完整的全端口信息来算！由于没有证明，也就无从考据这一点，只能对着实验结果说，可能是这样。

Balanced Truncation

https://web.stanford.edu/group/frg/course_work/CME345/CA-CME345-Ch6.pdf

假设电感矩阵可逆（虽然不太可能），重新将式（MOR）记为：

\begin{aligned} \dot{x_n}&=-\boldsymbol{C}^{-1}\boldsymbol{G}x_n+\boldsymbol{C}^{-1}\boldsymbol{B}u_p\\ &=\boldsymbol{\dot{A}}x_n+\boldsymbol{\dot{B}}u_p\\ i_p&=\boldsymbol{L^T}x_n \end{aligned} \tag{BTR a}

定义reachability (controllability) Gramian：

\mathcal{P}(t) = \int_{0}^{t} e^{\boldsymbol{\dot{A}} \tau} \boldsymbol{\dot{B}} \boldsymbol{\dot{B}}^T e^{\boldsymbol{\dot{A}}^* \tau} \, d\tau

定义observability Gramian（注意这里L原本就有转置）：

\mathcal{Q}(t) = \int_{0}^{t} e^{\boldsymbol{\dot{A}}^* \tau} \boldsymbol{L} \boldsymbol{L}^T e^{\boldsymbol{\dot{A}} \tau} \, d\tau

||x||_{\mathcal{A}}:=\sqrt{x^T\mathcal{A}x}称作semi-norm，可以用来衡量“能量”（为什么？），所以，如果能通过挑选合适的基，就可以消除符合以下定义的状态变量：

难以达到（difficult to reach）：这里指的是需要较大能量（||\mathbf{x}||_{\mathcal{P}^{-1}}^2）才能达到的状态。这意味着在系统中需要较大的输入才能使系统到达这些状态。
难以观测（difficult to observe）：指的是这些状态产生的观测能量较小（||\mathbf{x}||_{\mathcal{Q}}^2），即从输出信号中难以检测或观测到这些状态。

最终，我们希望找到一个基，使得状态的可达性和可观测性在某种程度上达到一个平衡点，既能较容易地达到这些状态，同时也能较容易地观测到它们。如果存在一个\tilde x=Tx，那么（BTR a）变形为：

\begin{aligned} \dot{\tilde x_n}&=T\boldsymbol{\dot{A}}T^{-1}\tilde x_n+T\boldsymbol{\dot{B}}u_p\\ i_p&=\boldsymbol{L^T}x_n \end{aligned} \tag{BTR b}

注意到exp(T^{-1}AT)=T^{-1}exp(A)T，则\exp(T\boldsymbol{\dot{A}}T^{-1}t)T\boldsymbol{\dot{B}}=T\exp(\boldsymbol{\dot{A}}t)B，那我们容易推得：

\tilde{\mathcal P}=T\mathcal PT^T\quad \tilde{\mathcal Q}=T^{-T}\mathcal Q T^{-1}

所以我们的目标就是找到相似变换阵T，使得P、Q获得尽可能好的性质，BTR的策略如下：

Cholesky分解得到\mathcal P = UU^T
特征值分解得到U^T\mathcal Q U=K\Sigma^2K^T
只需要令T=\Sigma^{1/2}K^{T}U^{-1}
有T\mathcal P T^T=T^{-T}\mathcal Q T^{-1}=\Sigma

误差分析（待看）

TICER

时域分析

https://quan-zhou.weebly.com/uploads/3/1/1/3/31137783/_3.pdf——系统的时域分析，不知道哪里的课件。

http://www.tup.tsinghua.edu.cn/upload/books/yz/077499-03.pdf——第 3 章线性系统的时域分析法，也不知道是清华出版社的哪本书。

假设一个星形电路，中心是节点N，N个端口标记为0到N-1。每一个分支由一个导电率 g_{iN}和一个电容c_{iN}组成，连接第i个端口和中心节点。当第i个端口施加一个阶跃电压，其他所有端口接地时，中心节点的KCL方程为：

\sum_{j\neq i}\left(g_{jN}v_N(t)+c_{jN}\frac{d v_N(t)}{dt}\right)=g_{iN}(u(t)-v_N(t))+c_{jN}\frac{d u(t)-v_N(t)}{dt}\\ \Rightarrow \dot v_N(t) + \frac{\gamma_N}{\chi_N}v_N(t) = \frac{g_{iN}}{\chi_N}u(t)+\frac{c_{iN}}{\chi_N}\delta (t)

记：

\gamma_N = \sum_{k=0}^{N-1} g_{kN}, \quad \chi_N = \sum_{k=0}^{N-1} c_{kN}, \quad \tau_N = \frac{\chi_N}{\gamma_N}

对冲激信号带入特解Ae^{-t/\tau_N}得到A={c_{iN}}/{\chi_N}；对阶跃信号设常数特解得到{g_{iN}}/{\gamma_N}，令响应y=Ce^{-t/\tau_N}+\frac{g_{iN}}{\gamma_N}，带入0初始状态得到C=-{g_{iN}}/{\gamma_N}，则总响应表示为：

h_{iN}(t) = \frac{g_{iN}}{\gamma_N} + \left(\frac{c_{iN}}{\chi_N} - \frac{g_{iN}}{\gamma_N}\right) e^{-t/\tau_N}

其他参考文档

pyMOR v2024.2.0.dev0 Manual

Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide. SIAM, Philadelphia, 2000.