Parameter Estimation for a Complex PDE System

Objective

This example will show you how to solve a complex system of PDEs, then optimize the system parameters to obtain excellent fit with experimental data. We will first obtain a solution using default values for the parameters with PDASOLVE (Figure a), then optimize the model using Excel Solver or NLSOLVE (Figure b).

Figure a
Figure b

3-Region coupled PDE System

For a detailed description of this problem, please see Example 1 of PDASOLVE

m11 Φ1t + m12 Φ2t = σx Φ1,xx - a x Φ1 - Φ2 - f t,x m21 Φ1t + m22 Φ2t + m23 Φ2t = κ x   Φ2,xx + a x Φ1 - Φ2 - f t,x m32 Φ2t + m33 Φ3t = σ x   Φ3,xx - a x ( Φ3 - Φ2 - f t,x

Properties

σx=10xlc0lc<x<ls1lsxla,      κx=β*10xlc0lc<x<ls1lsxla,      ax=γ1060xlc0lc<x<ls106lsxla,      ft,x=3e-t0xlc0else

Parameters

lc=0.0002,      β=1,      ls=0.0002254,      γ=0.5,      la=0.0004254,      mc=106,      mij=106-mc0-mc106-mc0-mc106

Initial Values

Φ1(0,x)= Φ2(0,x)= Φ3(0,x)=0

Boundary conditions

x = 0 x = lc x = ls x = la
σ0 Φ1,x=1 Φ1,x=0 Φ3,x=0 Φ3=0
Φ2,x=0 κlcΦ2,x- = κlc Φ2,x+ κls Φ2,x-= κls Φ2,x+ Φ2,x=0

Solution

Working with named system variables, parameters, and property functions we define the complete Excel model including equations, regions, and boundary conditions formulas as shown in Table 1.

Table 1
ABCDEFGH
1System Variables with initial valuesSystem Parameters
2tlc2.00E-04
3xls2.254E-04
4phi_10la4.254E-04
5phi_20beta1
6phi_30gamma0.5
7phi1xmc1E+06
8phi2x
9phi3xMass Matrix
10phi1xx1E+06=-mc0
11phi2xx=-mc1E+06=-mc
12phi3xx0=-mc1E+06
13
14Property Functionsx Domain0=la
15sigma=IF(x<=lc,1,IF(x<ls,0,1))t Interval01Boundary conditions matrix
16kappa=beta*IF(x<=lc,1,IF(x<ls,2,1))0R=sigma*phi1x-1
17a=gamma*IF(x<=lc,1E6,IF(x<ls,0,1E6))0N=phi2x
18f=IF(x<=lc,3*EXP(-2*t),0)=lcN=phi1x
19=lcC=kappa*phi2x
20System RHS EquationsRegions=lsC=kappa*phi2x
21eq1=sigma*phi1xx-a*(phi_1-phi_2-f)0=lc=lsN=phi3x
22eq2=kappa*phi2xx+a*(phi_1-phi_2-f)0=la=laN=phi2x
23eq3=sigma*phi3xx-a*(phi_3-phi_2-f)=ls=la=laD=phi_3

To aid readability, we assign assign names for selected Excel ranges in Table 1 which represent input arguments to PDASOLVE as shown in Table 2.

Table 2
Excel RangeAssigned Name
B21:B23Eqns
B2:B12Vars
F16:H23BCs
D14:E14xDom
D15:E15tDom
C10:E12M
C21:D23Rgns

Next, we evaluate the array formula =PDASOLVE(Eqns, Vars, BCs, xDom, tDom, M, Rgns, , {"FORMAT","TCOL1"}) in array G15:S27 and obtain the solution as shown in Table 3.

Table 3
GHIJKLMNOPQRS
15

x

0000.00020.00020.000200.00022540.00022540.00042540.00042540.00043
16

t

phi_1

phi_2

phi_3

phi_1

phi_2

phi_3

phi_1

phi_2

phi_3

phi_1

phi_2

phi_3

170000000000000
180.10.0490406-0.003126900.0491381-0.002664600-0.0026093-0.00039370-0.0022350
190.20.0888748-0.00478200.0889728-0.004394500-0.0043481-0.00033110-0.0040330
200.30.1204078-0.006139300.1205061-0.005812900-0.0057737-0.00028010-0.0055070
210.40.1451792-0.007253900.1452778-0.006977400-0.0069442-0.00023850-0.0067160
220.50.1644121-0.008169700.164511-0.007933900-0.0079055-0.00020460-0.007710
230.60.1791579-0.008924400.179257-0.008721700-0.0086973-0.00017690-0.0085280
240.70.1902228-0.009546800.190322-0.009371100-0.0093499-0.00015440-0.0092010
250.80.1983065-0.010062100.1984059-0.009908300-0.0098896-0.00013610-0.0097580
260.90.2039649-0.0104900.2040644-0.01035400-0.0103374-0.00012130-0.010220
2710.2076535-0.010846800.2077531-0.010725300-0.0107105-0.00010930-0.0106050

In Figure 1, we plot phi2(x,t) at the interface x = lc together with imperical values for Φ(lc,t). Clearly the default model does not agree well with experimental data. Next will show you how to optimize the model parameters so it agrees with experimental data accurately.

Figure 1

Parameter Estimation

The observed experimental values for Φ(lc,t) are give in Table 6.

Table 6
MN
1timeEmperical phi2
200
30.1-0.002052705
40.2-0.00380538
50.3-0.005944579
60.4-0.006668727
70.5-0.008436437
80.6-0.009136398
90.7-0.009734491
100.8-0.010031262
110.9-0.010050529
121-0.010387955

Our task is to compute optimal values for the design parameters γ, and mc which minimize the sum of square residuals between our model predictions and the observed data in Table 6. We have two options to solve this dynamical minimization problem: Excel native Solver, and NLSOLVE. We demonstrate both procedures below.

Using Excel Solver

We define an objective formula to be minimized. The objective formula calculates the summation of the square residuals between the observed values and their corresponding computed values from the solution result in Table 3. The definition of the objective formula is shown in S2 of Table 7 with initial value of 0.00033.

Table 7
PS
1=(L18-N2)^2Objective
2=(L19-N3)^2=SUM(P1:P10)
3=(L20-N4)^2
4=(L21-N5)^2
5=(L22-N6)^2
6=(L23-N7)^2
7=(L24-N8)^2
8=(L25-N9)^2
9=(L26-N10)^2
10=(L27-N11)^2
10=(L28-N12)^2

The remaining task is to configure and run Excel Solver command. We configure the solver exactly as shown in Figure 2. We select to minimize the objective formula S2, by varying gamma and mc, subject to the upper bound constraint mc <= 10e6 which is added for numerical stability since larger values for the off-diagonal mass matrix element mc leads to a divergent system.

Figure 2

We run the solver which reports a feasible solution in less than a minute of computing time as shown in the Answer Report generated by the Excel Solver in Figure 3. Upon accepting the solution, Excel automatically updates the spreadsheet to reflect the optimal computed values. The optimized solution for phi2 is shown in Figure 4 exhibiting excellent agreement with observed values.

Figure 3
Figure 4

Using NLSOLVE

Unlike Excel Solver which requires an objective formula, NLSOLVE requires the list of individual residual constraint formulas which we define as shown in Table 8.

Table 8
R
1=ARRAYVAL(L18)-N2
2=ARRAYVAL(L19)-N3
3=ARRAYVAL(L20)-N4
4=ARRAYVAL(L21)-N5
5=ARRAYVAL(L22)-N6
6=ARRAYVAL(L23)-N7
7=ARRAYVAL(L24)-N8
8=ARRAYVAL(L25)-N9
9=ARRAYVAL(L26)-N10
10=ARRAYVAL(L27)-N11
10=ARRAYVAL(L27)-N12

What is ARRAYVAL? ARRAYVAL is one of the criterion functions that allow us to define dynamic constraints on the solution array. In this case we are simply picking a value from the solution but in general, a criterion function can be used to specify more complex constraint on the solution.

Next, we evaluate the formula =NLSOLVE(R1:R10, (gamma, mc)) in array S4:T8 and obtain virtually the same solution found by Excel Solver as shown in Table 9. However, NLSOLVE, converges in fewer iterations in less than one third of the time. Since NLSOLVE is a pure spreadsheet function, it does not modify its arguments. Therefore, to update the spreadsheet result, the computed values for gamma and mc are copied manually to their respective variable cells D6 and D7 of Table 3 of Example 1.

Table 9
ST
4gamma0.195638747
5mc958050.8287
6SSERROR1.61147E-06
7ITRN10
8TIME (s)18.729

References

Ghaddar, C.K. "Rapid Modeling and Parameter Estimation of Partial Differential Algebraic Equations by a Functional Spreadsheet Paradigm". Math. Comput. Appl. 2018, 23, 39.
Available at: http://www.mdpi.com/2297-8747/23/3/39

Question or Comment? Email us:
support @ excel-works.com
ExceLab: Transforming Excel into a Calculus Power House
ExceLab functions and methods are protected by USA Patents 9286286, 9892108, 10114812 and pending.
© 2015-2019, ExcelWorks LLC
Boston, USA