WIP: Prototyping Pyomo.DoE lite using symbolic differentiation to assemble sensitivity matrix

[adowling2]

This is inspired by Pyomo/pyomo#3564

[adowling2]

@michaelbynum Thank you for your help.

@djlaky Check this out. I got a proof-of-concept "lite" Pyomo.DoE using that uses symbolic differentiation instead of finite difference. See skunk_works/doe_lite.ipynb. This needs a lot of work to be generalized, but it runs great on the one test problem with A-optimality. D-optimality using Cholesky factorization should not be that bad.

[adowling2]

I successfully prototyped A- and D-optimality on a toy problem. The next step is to try it with the TCLab.

[adowling2]

Updates:

• Moved "DoE lite" to a .py file
• Confirmed same optimal FIM as original Pyomo.DoE for the reactor problem
• Started to debug DoE lite with TCLab:
  • Restoration failure for D-optimality
  • Lots of line search for A-optimality

Here are the first 10 iterations for A-optimality:

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -3.5558612e+08 3.20e-04 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -3.7981769e+08 4.28e+05 1.77e+00  -1.0 7.79e+12    -  2.31e-05 2.86e-06f  3
   2 -3.8336197e+08 4.39e+05 1.76e+00  -1.0 8.54e+12    -  3.40e-07 4.24e-07f  7
   3 -3.8688642e+08 4.50e+05 1.72e+00  -1.0 8.78e+12    -  1.23e-08 4.13e-07f  7
   4 -3.9038854e+08 4.60e+05 1.71e+00  -1.0 8.98e+12    -  2.93e-07 4.03e-07f  7
   5 -3.9125525e+08 4.61e+05 1.71e+00  -1.0 9.22e+12    -  4.31e-08 9.75e-08f  9
   6 -3.9211961e+08 4.62e+05 1.71e+00  -1.0 9.28e+12    -  1.28e-07 9.67e-08f  9
   7 -3.9233466e+08 4.62e+05 1.72e+00  -1.0 9.41e+12    -  8.23e-08 2.38e-08f 11
   8 -3.9244210e+08 4.62e+05 1.72e+00  -1.0 9.44e+12    -  7.89e-08 1.19e-08f 12
   9 -3.9254956e+08 4.62e+05 1.73e+00  -1.0 9.43e+12    -  1.00e-07 1.19e-08f 12
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -3.9260332e+08 4.62e+05 1.80e+00  -1.0 9.42e+12    -  8.42e-07 5.96e-09f 13

Number of Iterations....: 10

Here are some key results from the diagnostics toolbox:

The following variable(s) have extreme values (<1.0E-04 or > 1.0E+04):

    Ts1dot[456.02]: 3.453737580852642e-05
    Ts1dot[600.01]: 1.2155806483477504e-06
    jac_variables_wrt_param[1,605]: 8.976067993358789e-05
    jac_variables_wrt_param[76,605]: -1.3898542405420364e-05
    jac_variables_wrt_param[151,605]: -9.9963016122793e-05
    jac_variables_wrt_param[154,605]: 1.5515115174230433e-05
    jac_variables_wrt_param[220,605]: 7.215000065399986e-05
    jac_variables_wrt_param[223,605]: -3.0306219981987373e-05
    jac_variables_wrt_param[301,605]: -8.968663966169215e-06
    jac_variables_wrt_param[332,605]: 5.264566509083845e-05
    jac_variables_wrt_param[373,605]: 4.515108544436486e-05
    jac_variables_wrt_param[376,605]: -6.62738958357418e-05
    jac_variables_wrt_param[448,605]: -9.856574061871153e-05
    fim[602,602]: 332820317.93551016
    fim[602,603]: 4329226.052524379
    fim[602,604]: -10061223.441862779
    fim[603,602]: 4329226.052524379
    fim[603,603]: 57277204.12901004
    fim[603,604]: 8297477.087120233
    fim[603,605]: 97316.47877611648
    fim[604,602]: -10061223.441862779
    fim[604,603]: 8297477.087120234
    fim[604,604]: 2505630.9302714965
    fim[604,605]: 15302.763696236507
    fim[605,603]: 97316.47877611648
    fim[605,604]: 15302.763696236505

The following constraint(s) are associated with extreme Jacobian values (<1.0E-04 or>1.0E+04):

    fim_constraint[602,602]: 7.302E+04
    fim_constraint[602,603]: 3.953E+04
    fim_constraint[603,602]: 3.953E+04
    fim_constraint[602,604]: 3.665E+04
    fim_constraint[604,602]: 3.665E+04
    fim_constraint[602,605]: 3.651E+04
    fim_constraint[605,602]: 3.651E+04
    fim_constraint[603,603]: 3.029E+04
    fim_constraint[603,604]: 1.547E+04
    fim_constraint[604,603]: 1.547E+04
    fim_constraint[603,605]: 1.514E+04
    fim_constraint[605,603]: 1.514E+04

The following constraint(s) have large residuals (>1.0E-05):

    Th1_ode[216.01]: 5.30907E-05
    Th1_ode[222.0]: 5.30907E-05
    Th1_ode[228.01]: 5.30907E-05
    Th1_ode[234.01]: 5.30907E-05
    Th1_ode[516.02]: 5.30907E-05
    Th1_ode[522.01]: 5.30907E-05
    Th1_ode[528.0]: 5.30907E-05
    Th1_ode[534.01]: 5.30907E-05
    Th1_ode[816.01]: 5.30907E-05
    Th1_ode[822.01]: 5.30907E-05
    Th1_ode[828.0]: 5.30907E-05
    Th1_ode[834.0]: 5.30907E-05
    jacobian_constraint[36,604]: 3.19998E-04
    jacobian_constraint[37,604]: 3.19998E-04
    jacobian_constraint[38,604]: 3.19998E-04
    jacobian_constraint[39,604]: 3.19998E-04
    jacobian_constraint[86,604]: 3.19998E-04
    jacobian_constraint[87,604]: 3.19998E-04
    jacobian_constraint[88,604]: 3.19998E-04
    jacobian_constraint[89,604]: 3.19998E-04
    jacobian_constraint[136,604]: 3.19998E-04
    jacobian_constraint[137,604]: 3.19998E-04
    jacobian_constraint[138,604]: 3.19998E-04
    jacobian_constraint[139,604]: 3.19998E-04
    fim_constraint[602,602]: 4.61733E+05
    fim_constraint[602,603]: 2.17456E+04
    fim_constraint[602,604]: 7.29394E+03
    fim_constraint[602,605]: 2.72927E+01
    fim_constraint[603,602]: 2.17456E+04
    fim_constraint[603,603]: 6.18722E+04
    fim_constraint[603,604]: 1.08222E+04
    fim_constraint[603,605]: 1.06842E+02
    fim_constraint[604,602]: 7.29394E+03
    fim_constraint[604,603]: 1.08222E+04
    fim_constraint[604,604]: 3.46261E+03
    fim_constraint[604,605]: 2.01390E+01
    fim_constraint[605,602]: 2.72927E+01
    fim_constraint[605,603]: 1.06842E+02
    fim_constraint[605,604]: 2.01390E+01
    fim_constraint[605,605]: 1.85838E-01

Ideas:

• Try scaling the FIM (number of measurements?)
• Try Degeneracy Hunter
• Try SVD toolbox

[adowling2]

More ideas:

• Add scaling
• Check Jacobian against Pyomo.DoE