INVERSION: CALIBRATION OF MULTIPLE FORWARD MODELS¶

This example shows how multiple computational models can be calibrated simultaneously, a process known in some fields as joint inversion. In this example, two different sets of deformation measurements of the same specimen are used to calibrate its Young's modulus $E$, by subjecting it to an uncertain distributed load $p$ and to an uncertain point load $P$.

The first set of measurements refers to the mid-span deflection of the simply supported specimen under the distributed load $p$, while and the second set refers to its elongation under the constant load $P$.

Package imports¶

In [1]:
from uqpylab import sessions
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image

Start a remote UQCloud session¶

In [2]:
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

Processing .
.
 done!

 uqpylab.sessions :: INFO     :: This is UQ[py]Lab, version 1.00, running on https://uqcloud.ethz.ch. 
                                 UQ[py]Lab is free software, published under the open source BSD 3-clause license.
                                 To request special permissions, please contact:
                                  - Stefano Marelli (marelli@ibk.baug.ethz.ch).
                                 A new session (ec38a1d451564d789f01d03ae1f48511) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [3]:
uq.rng(100, 'twister');

Forward models¶

The first forward model is the deflection of the beam at the mid-span location ($V$), which reads:

$$ V = \frac{ 5 p L^4 }{32 E b h^3}$$

This computation is carried out by the function SimplySupportedBeam(X). The input variables of this function are gathered into the $N \times M$ matrix X, where $N$ and $M$ are the number of realizations input variables, respectively. The variables are given in the following order:

  • $b$: beam width $(m)$
  • $h$: beam height $(m)$
  • $L$: beam length $(m)$
  • $E$: Young's modulus $(MPa)$
  • $p$: uniform load $(kN/m)$

The simply supported beam problem is shown in the following figure:

In [4]:
Image(filename='SimplySupportedBeam.png',width=800)
Out[4]:
No description has been provided for this image

Create a UQ[py]Lab MODEL object:

In [5]:
ModelOpts1 = {
    'Type' : 'Model',
    'Name': 'Beam mid-span deflection',
    'mString': '5/32.*X(:,5).*X(:,3).^4./(X(:,4).*X(:,1).*X(:,2).^3)',
    'isVectorized': 1
}
In [6]:
myModel1 = uq.createModel(ModelOpts1)

The second forward model is the elongation of a prismatic beam ($U$), which is computed by:

$$ U = \frac{PL}{Ebh}$$

with an additional parameter:

  • $P$: Point load $(N)$

The beam elongation problem is shown in the following figure:

In [7]:
Image(filename='ElongationBeam.png',width=800)
Out[7]:
No description has been provided for this image

Create a UQ[py]Lab MODEL object.

The model in this case is directly defined with a string:

In [8]:
ModelOpts2 = {'Type': 'Model', 
             'Name': 'Beam elongation',
             'mString':'X(:,5).*X(:,3)./(X(:,1).*X(:,2).*X(:,4))',
             'isVectorized': 1}
In [9]:
myModel2 = uq.createModel(ModelOpts2)

For both forward models, store them in a list of dictionaries.

In addition, specify a |PMap| vector to specify which prior model parameters are to be passed to the first forward model.

The second model shares most inputs with the previous, but instead of the uniform load $p$ (parameter |5|), it uses the point load $P$ (parameter |6|). This is specified with the |PMap| vector in the second dictionary:

In [10]:
forwardModels = [
    {
        'Model': myModel1['Name'],
        'PMap': [1,2,3,4,5]
    },
    {
        'Model': myModel2['Name'],
        'PMap': [1,2,3,4,6]
    }
]

Prior distribution of the model parameters¶

The prior distribution of the model parameter $E$ and the uncertain loads $p$ and $P$ are given by:

  • $E \sim \mathcal{LN}(\mu_E = 30\times10^9, \sigma_E = 4.5\times10^9)~(N/m^2)$
  • $p \sim \mathcal{N}(\mu_p = 12\times10^3, \sigma_p = 3\times10^3)~(N/m)$
  • $P \sim \mathcal{LN}(\mu_P = 50\times10^3, \sigma_P = 10^3)~(N)$

The parameters $b$, $h$, and $L$, on the other hand, are fixed as follows:

  • $b = 0.15~(m)$
  • $h = 0.3~(m)$
  • $L = 5~(m)$

Gather these distributions in an INPUT object:

In [11]:
PriorOpts = {
  "Marginals": [
    {
        "Name": "b", # beam width
        "Type": "Constant",
        "Parameters": [0.15] # (m)
    },
    {
        "Name": "h", # beam height
        "Type": "Constant",
        "Parameters": [0.3] # (m)
    },
    {
        "Name": "L", # beam length
        "Type": "Constant",
        "Parameters": [5] # (m) 
    },
    {
        "Name": "E", # Young's modulus
        "Type": "LogNormal",
        "Moments": [30e9,4.5e9] # (N/m^2)
    },
    {
        "Name": "p", # uniform load
        "Type": "LogNormal",
        "Moments": [12000,3000] # (N/m)
    },
    {
        "Name": "P", # point load
        "Type": "LogNormal",
        "Moments": [50e3, 1e3] # (N)
    }
  ]
}
In [12]:
myPriorDist = uq.createInput(PriorOpts)

Measurement data¶

In the case of multiple forward models, and therefore different types of data, it is necessary to define the full |MOMap| list to identify which model output needs to be compared with which data set:

In [13]:
V_mid = np.array([12.84, 13.12, 12.13, 12.19, 12.67])/1000 # (m)
U_right = np.array([0.235, 0.236, 0.229])/1000             # (m)

myData = [
    # Data group 1
    {
        'y': V_mid.tolist(),
        'Name': 'Beam mid-span deflection',
        'MOMap': [1,        # Model ID
                  1]        # Output ID
    },
    # Data group 2
    {
        'y': U_right.tolist(),
        'Name': 'Beam elongation',
        'MOMap': [2,        # Model ID
                  1  ]      # Output ID
    }
]

Discrepancy model¶

To infer the discrepancy variance, uniform priors are put on the discrepancy parameters:

  • $\sigma^2_{\mathrm{1}} \sim \mathcal{U}(0, 10^{-6})~(m^2)$
  • $\sigma^2_{\mathrm{2}} \sim \mathcal{U}(0, 10^{-10})~(m^2)$

Create two UQ[py]Lab INPUT objects based on these distributions:

In [14]:
SigmaOpts = {
    'Marginals': [
        {
            'Name': 'Sigma2V',
            'Type': 'Uniform',
            'Parameters': [0,1e-6]      # (m^2)
        }
    ]
}

SigmaDist1 = uq.createInput(SigmaOpts)

SigmaOpts = {
    'Marginals': [
        {
            'Name': 'Sigma2U',
            'Type': 'Uniform',
            'Parameters': [0,1e-10]    # (m^2)
        }
    ]
}

SigmaDist2 = uq.createInput(SigmaOpts)

Assign these distributions to the discrepancy model options:

In [15]:
DiscrepancyOpts = [
    {
        'Type': 'Gaussian',
        'Prior': SigmaDist1['Name']
    },
    {
        'Type': 'Gaussian',
        'Prior': SigmaDist2['Name']
    }
]

Bayesian analysis¶

MCMC solver options¶

To sample from the posterior distribution, the affine invariant ensemble algorithm is employed in this example, with $200$ iterations and $100$ parallel chains:

In [16]:
Solver = {
    'Type': 'MCMC',
    'MCMC': {
        'Sampler': 'AIES',
        'Steps': 200,
        'NChains': 100
    }
}

6.2 Posterior sample generation¶

The options of the Bayesian analysis are specified with the following dictionary, where the |forwardModels| list of dictionaries is passed to the |BayesOpts['ForwardModel']| key-value pair:

In [17]:
BayesOpts = {
    "Type": "Inversion",
    "Name": "Bayesian model1",
    "Prior": myPriorDist['Name'],
    "ForwardModel": forwardModels,
    "Data": myData,
    "Discrepancy": DiscrepancyOpts,
    "Solver": Solver
}

Run the Bayesian inversion analysis:

In [18]:
myBayesianAnalysis = uq.createAnalysis(BayesOpts)

Processing .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
 done!

Print out a report of the results:

In [19]:
uq.print(myBayesianAnalysis)

%----------------------- Inversion output -----------------------%
   Number of calibrated model parameters:         3
   Number of non-calibrated model parameters:     3

   Number of calibrated discrepancy parameters:   2

%------------------- Data and Discrepancy
%  Data-/Discrepancy group 1:
   Number of independent observations:            5

   Discrepancy:
      Type:                                       Gaussian
      Discrepancy family:                         Scalar
      Discrepancy parameters known:               No

   Associated outputs:
      Model 1: 
         Output dimensions:                       1

%  Data-/Discrepancy group 2:
   Number of independent observations:            3

   Discrepancy:
      Type:                                       Gaussian
      Discrepancy family:                         Scalar
      Discrepancy parameters known:               No

   Associated outputs:
      Model 2: 
         Output dimensions:                       1

%------------------- Solver
   Solution method:                               MCMC

   Algorithm:                                     AIES
   Duration (HH:MM:SS):                           00:01:15
   Number of sample points:                       2.00e+04

%------------------- Posterior Marginals
---------------------------------------------------------------------
| Parameter | Mean    | Std     | (0.025-0.97) Quant. | Type        |
---------------------------------------------------------------------
| E         | 2.4e+10 | 6.5e+08 | (2.3e+10 - 2.5e+10) | Model       |
| p         | 1.3e+04 | 4.4e+02 | (1.2e+04 - 1.3e+04) | Model       |
| P         | 5e+04   | 1.1e+03 | (4.8e+04 - 5.2e+04) | Model       |
| Sigma2V   | 4.3e-07 | 2.5e-07 | (9.2e-08 - 9.5e-07) | Discrepancy |
| Sigma2U   | 4.5e-11 | 2.6e-11 | (7.1e-12 - 9.6e-11) | Discrepancy |
---------------------------------------------------------------------

%------------------- Point estimate
----------------------------------------
| Parameter | Mean    | Parameter Type |
----------------------------------------
| E         | 2.4e+10 | Model          |
| p         | 1.3e+04 | Model          |
| P         | 5e+04   | Model          |
| Sigma2V   | 4.3e-07 | Discrepancy    |
| Sigma2U   | 4.5e-11 | Discrepancy    |
----------------------------------------

%------------------- Correlation matrix (model parameters)
------------------------------
|   |  E       p       P     |
------------------------------
| E |  1       0.78    0.79  |
| p |  0.78    1       0.61  |
| P |  0.79    0.61    1     |
------------------------------

%------------------- Correlation matrix (discrepancy parameters)
----------------------------------
|         |  Sigma2V    Sigma2U  |
----------------------------------
| Sigma2V |  1          -0.043   |
| Sigma2U |  -0.043     1        |
----------------------------------

Create a graphical representation of the results:

In [20]:
uq.display(myBayesianAnalysis);

Terminate the remote UQCloud session¶

In [21]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session ec38a1d451564d789f01d03ae1f48511 terminated.
Out[21]:
True