天天看點

PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)

💥💥💥💞💞💞歡迎來到本部落格❤️❤️❤️💥💥💥

📋📋📋本文目錄如下:⛳️⛳️⛳️

目錄

​​1 概述​​

​​2 數學模型​​

​​3 算例及仿真結果​​

​​4 Matlab代碼及文章詳細閱讀​​

​​5 結論​​

​​6 寫在最後​​

1 概述

在逆問題中,目标是從間接觀察中推斷實體參數(例如密度、聲速或電導率)。當基礎模型由偏微分方程 (PDE)(例如,波動方程或麥克斯韋方程)描述時,觀察到的資料通常是多個右手邊 PDE 解的部分測量值。這些參數通常顯示為 PDE 中的系數。這些問題出現在許多應用中,例如地球實體學 [1、2、3、4]、醫學成像 [5、6] 和無損檢測。

許多逆和參數估計問題可以寫成偏微分方程限制的優化問題。是以,目标是從幾個右側偏微分方程解的部分測量中推斷出參數,通常是偏微分方程的系數。這種偏微分方程限制的問題可以通過找到拉格朗日量的靜止點來解決,這需要同時更新參數和(伴随)狀态變量。對于大規模問題,這種一次性全部的方法不可行,因為它需要存儲所有狀态變量。在這種情況下,人們通常采用簡化方法,其中通過求解偏微分方程來顯式消除限制(在每次疊代中)。這兩種方法及其變體是解決由逆問題引起的偏微分方程限制優化問題的主要主力。在本文中,我們提出了一種替代方法,旨在結合這兩種方法的優點。我們的方法基于限制優化問題的二次懲罰公式。通過消除狀态變量,我們開發了一種有效的算法,該算法具有與傳統簡化方法大緻相同的計算複雜性,同時利用更大的搜尋空間。數值結果表明,該方法确實降低了問題的一些非線性,并且對初始疊代的敏感性較低。

2 數學模型

對于線性偏微分方程,可以将逆問題(離散化後)表述為以下形式的限制優化問題:

PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)

在實踐中,通常在公式 (1) 中包含一個正則化項,以減輕問題的不适定性。然而,為了簡化讨論,忽略了這些術語,因為可以在需要時添加适當的正則化術語。

在由逆問題引起的應用中,限制問題 (1) 通常使用拉格朗日乘數 [7, 8] 或順序二次規劃 (SQP) [9, 10] 的方法來解決。這需要同時優化參數、狀态和拉格朗日乘數(或伴随狀态變量)。雖然從優化的角度來看,這種一次性方法通常非常有吸引力,但它們通常不适用于大規模問題,因為我們無法同時存儲所有 K 個實驗的狀态變量。相反,所謂的簡化方法是基于(塊)消除限制來制定參數上的無限制優化問題:

PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)

如果對狀态進行完整測量(即 P 是可逆的),則所述逆問題将更容易解決。在這種情況下,我們可以從資料中重建狀态為 u = P−1d,然後通過求解來恢複參數:

PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)

詳細數學模型見第4部分。

3 算例及仿真結果

PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)
PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)
PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)
PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)
PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)
PDE優化|逆問題中偏微分方程限制優化的懲罰方法(Matlab代碼實作)

 4 Matlab代碼及文章詳細閱讀

n  = 51;
f  = 5;
x  = linspace(0,1,n-1);
h  = x(2) - x(1);
mt = ones(n-1,1);
omega = 2*pi*f;

D  = spdiags(ones(n,1)*[-1 1]/h,[0:1],n-1,n);
w  = ones(n,1); w([1 end]) = .5;
At = (1i*omega)*diags(w) + D'*diags(mt)*D;

[Vr,Dr] = eig(full(At'*At)); [Dr,Ir] = sort(diag(Dr),'descend'); Vr = Vr(:,Ir);
kappar  = Dr(1)/Dr(end);

a  = [1e-1 1 1e1];
L  = [1 10 20];
ls = {'ro','bx','g^'};

P{1} = speye(n);
P{2} = fft(eye(n))/sqrt(n);

for k = 1:1
    for l = 1:length(L)  
        %I = randperm(n); I  = I(1:L(l));
        I = round(linspace(1,n,L(l)));
        Pl = P{k}(:,I);

        mu = eigmax(@(x)At'\(Pl*Pl'*(At\x)),n);

        figure;semilogy(1:n,Dr,'k*','markersize',10);hold on;
        xlabel('n');ylabel('特征值');xlim([1 n]);

        for i = 1:length(a)
            lambda = a(i)*real(mu);
            B = (At'*At) + (1/lambda)*(Pl*Pl');
            %[Vp,Dp] = eig(full(diags(1./Dr)*ex(Vr'*B*Vr))); [Dp,Ip] = sort(diag(Dp),'descend'); Vp = Vp(:,Ip);
            [Vp,Dp] = eig(full(B)); [Dp,Ip] = sort(diag(Dp),'descend'); Vp = Vp(:,Ip);
            kappap(l,i) = Dp(1)/Dp(end);
            lower = 1/(1+L(l)/(lambda*Dr(end)));
            upper =   (1+L(l)/(lambda*Dr(1)));
            
            semilogy(1:n,abs(Dp),ls{i},'markersize',10);hold on;
        end
        legend('reduced',['\lambda = 0.1'],['\lambda = 1.0'],['\lambda = 10'],'location','SouthWest');
    end

    kappa{k} = kappap/kappar;
end

savefig(1:3,'../../doc/figs/example2');

latextable(kappa{1},'Horiz',{'$\lambda = 0.1$','$\lambda = 1$','$\lambda = 10$'},'Vert',{['L = ' num2str(L(1))],['L = ' num2str(L(2))],['L = ' num2str(L(3))] },'Hline',[1 NaN],'format','%1.2e','name','../../doc/figs/example2.tex');

      

 5 結論

6 寫在最後

繼續閱讀