Fitting data with interdependence

Question

I have an issue, of which I'm not sure it can be solved, with or without Mathematica. I feel as if it should be possible, but I'm quite clueless as to how.

Some introduction is required. I'm working in a research lab, and we are trying to characterize some samples that were fabricated. These are chips with capacitances on them, three different types. Lets call them $C_1$, $C_2$ and $C_J$. What we then have is 5 devices on which $C_1$ and $C_2$ are constant and $C_J$ is varied, and 4 devices on which $C_J$ is constant and $C_1$ and $C_2$ are varied. However, this is fabrication, we are not sure of the exact values of these parameters. This is what the experiment wants to characterize. They are both being varied by a parameter which we do know: the number of fingers we write on the chip, lets call them $F_{1,2}$ and $F_J$. So the idea is to find the functions $C_{1,2}(F_{1,2})$ and $C_J(F_J)$. Moreover, and this is a little messy, in these measurements it is the case that $C_1$ and $C_2$ always consist of the same amount of fingers, they simply have a slightly different length and thus lead to different capacitance values.

On these devices, measurements are performed. However, we can't measure the capacitance directly; we measure some energy related parameters, $\omega_{1,2}$ and $J$. This is where it gets a bit tricky, because they have a complicated definition, which might be easiest to just give in code:

Cmatrix = ( {
    {C1 + CJ + CJ, -CJ, 0, -CJ},
    {-CJ, C2 + CJ + CJ, -CJ, 0},
    {0, -CJ, C1 + CJ + CJ, -CJ},
    {-CJ, 0, -CJ, C2 + CJ + CJ}
   } );

InvCmat = Inverse[Cmatrix];

ω1 = Sqrt[1/L InvCmat[[1, 1]]];
ω2 = Sqrt[1/L InvCmat[[2, 2]]];

J = 1/2 InvCmat[[1, 2]]/Sqrt[InvCmat[[1, 1]] InvCmat[[2, 2]]] Sqrt[ω1*ω2];

I should add that in the above L is a known parameter, with value 1.84.

So, what I have is essentially nine datapoints:

{{20, 20, 2, 7.70438, 7.50004, 0.156103}, {20, 20, 3, 7.60618, 
  7.30572, 0.212588}, {20, 20, 4, 7.46063, 7.18521, 0.273479}, {20, 
  20, 6, 7.25941, 6.79769, 0.354176}, {20, 20, 10, 6.93511, 6.59471, 
  0.494317}, {13, 13, 3, 8.8022, 8.20026, 0.297579}, {16, 16, 3, 
  8.22736, 7.78623, 0.252157}, {20, 20, 3, 7.58361, 7.359, 
  0.212081}, {24, 24, 3, 7.08458, 6.90351, 0.178109}}

Each each containing (in order of occurance) $F_1$, $F_2$, $F_J$, $\omega_1$, $\omega_2$, $J$. From these datapoints I'd like to find a model for $C_i(F_i)$ and $C_J(F_J)$. Maybe an exact model is not possible, but an approximation should be. Based on physical intuition the relationships should be linear, but it might get a little messy due to fabricational issues.

So, given all this, would anyone have an idea on how to start?

I should also add that the fact that $C_1 \neq C_2$ is an error in the fabrication. Ideally this would not be the case, which simplifies the problem for sure, but even then I'd not know how to do such a fitting routine. So if instead you want to set $C_1 = C_2$ in the initial matrix and describe a method from there it would also be valuable.

Edit: I forgot to add. All values of $C$ are positive and real and $probably$ on the order of $10^{-4}$.

Answer 1

You pose a two part question: Part 1 find the "C's" from omega and J:

Cmatrix = ({{C1 + CJ + CJ, -CJ, 0, -CJ}, {-CJ, C2 + CJ + CJ, -CJ, 
 0}, {0, -CJ, C1 + CJ + CJ, -CJ}, {-CJ, 0, -CJ, C2 + CJ + CJ}});
InvCmat = Inverse[Cmatrix];
L = 1.84
\[Omega]1 = Sqrt[1/L InvCmat[[1, 1]]];
\[Omega]2 = Sqrt[1/L InvCmat[[2, 2]]];
J = 1/2 InvCmat[[1, 2]]/
   Sqrt[InvCmat[[1, 1]] InvCmat[[2, 2]]] Sqrt[\[Omega]1*\[Omega]2];
 NMinimize[{({\[Omega]1, \[Omega]2, J}
      - {7.70438, 7.50004, 0.156103})^2 // Total, C1 > 0, C2 > 0, 
     CJ > 0}, {C1, C2, CJ}]

{7.89631*10^-31, {C1 -> 0.00841431, C2 -> 0.00892173, CJ -> 0.000386302}}

Here are plots of Ci vs Fi , (i=1,2,J)

Answer 2

One can use FindMinimum and an objective function to find parameters that minimize the fit between the data and reconstructed data (i.e., forward model) provided one has equations for the forward model.

This will be a bit long winded but here goes.

Step 1. Define the model

We will assume that the capacitor is a linear function of the finger data (you are free to use almost any model).

C1 = c1*F1 + b1

where C1 is the unknown capacitance, F1 is the known finger, and c1 and b1 are the slope and intercept.

We will use a matrix that incorporates this model. We use ReplaceAll to convert your original matrix into the model matrix.

cmodel = ({
    {C1 + CJ + CJ, -CJ, 0, -CJ},
    {-CJ, C2 + CJ + CJ, -CJ, 0},
    {0, -CJ, C1 + CJ + CJ, -CJ},
    {-CJ, 0, -CJ, C2 + CJ + CJ}
    }) /.
   {C1 -> c1*F1 + b1,
    C2 -> c2*F2 + b2,
    CJ -> cj*Fj + bj}

Which results in

{
  {b1 + c1 F1 + 2 (bj + cj Fj), -bj - cj Fj, 
   0, -bj - cj Fj}, {-bj - cj Fj, 
   b2 + c2 F2 + 2 (bj + cj Fj), -bj - cj Fj, 0}, {0, -bj - cj Fj, 
   b1 + c1 F1 + 2 (bj + cj Fj), -bj - cj Fj}, {-bj - cj Fj, 
   0, -bj - cj Fj, b2 + c2 F2 + 2 (bj + cj Fj)}
  }

Step 2. Define the forward model functions

invCmodel = Inverse[cmodel];

ω1 = Sqrt[1/1.84 invCmodel[[1, 1]]];

ω1rec[c1_, b1_, c2_, b2_, cj_, bj_, F1_, F2_, Fj_] = ω1;

ω2 = Sqrt[1/1.84 invCmodel[[2, 2]]];

ω2rec[c1_, b1_, c2_, b2_, cj_, bj_, F1_, F2_, Fj_] = ω2;

J = 1/2 invCmodel[[1, 2]] Sqrt[ω1*ω2]/
    Sqrt[invCmodel[[1, 1]] invCmodel[[2,2]]];

Jrec[c1_, b1_, c2_, b2_, cj_, bj_, F1_, F2_, Fj_] = J;

Note the use of Set (i.e, =) in defining the functions rather than the normal SetDelayed (i.e., :=).

Step 3. Define an objective function

We will take the measured data and compute reconstructed data. It is easy to make this more terse but a longer version may be easier to understand.

data is the set of measured data with F1 in column 1, ... ω1 in column 4 ...

The variable parameterVector represents the list {c1, b1, c2, b2, cj, bj}.

In the code funny lines that use the Sequence function enable us to use lists as arguments to the functions (e.g., ω1rec).

Reconstructed data is computed, then the residual is equated to the difference between the measured and reconstructed data. The objective function is the sum of the squares of the residuals.

objectiveFunction[data_, parameterVectors_] := Module[{
   dataReconstructed,
   residual
   },

  (* Apply the forward model equations to the
      finger data to get reconstructed measurements *)

  dataReconstructed = Map[Function[fingers,
     {
      ω1rec[Sequence @@ parameterVectors, Sequence @@ fingers],
      ω2rec[Sequence @@ parameterVectors, Sequence @@ fingers],
        Jrec[Sequence @@ parameterVectors, Sequence @@ fingers]
      }],
    data[[All, 1 ;; 3]]
    ];

  (* Compute the residual *)

  residual = dataReconstructed - data[[All, 4 ;; 6]];

  (* The objective function is the sum of
      the residuals squared *)

  Total@Total@Map[#^2 &, residual, {2}]

  ]

Step 4. Apply FindMinimum

I did this twice (I'll only show you the second time). I discovered that the offset for Cj was less than zero.

FindMinimum can use constraints and works best if you supply it with reasonable initial values.

The solution, sol, represents the objective function and values for the parameters provided as rules.

sol = FindMinimum[
  {
   objectiveFunction[data, {c1, b1, c2, b2, cj, bj}],
   c1 > 0 && c1 ∈ Reals && b1 > 0 && b1 ∈ Reals &&
    c2 > 0 && c2 ∈ Reals && b2 > 0 && b2 ∈ Reals &&  
    cj > 0 && cj ∈ Reals && bj > -1*10^-4 && bj ∈ Reals
   },
  {{c1, 3.4*10^-4}, {b1, 1.6*10^-3}, {c2, 3*10^-4}, {b2, 
    3*10^-3}, {cj, 1.9*10^-4}, {bj, -4*10^-5}}
  ]

the result is

{0.0539887, {c1 -> 0.000340998, b1 -> 0.00157195, c2 -> 0.000301931, 
  b2 -> 0.00316187, cj -> 0.000192957, bj -> -0.0000434626}}

Step 5. Check the results

We will generate the reconstructed data to see how it compares to the measured data. I elected to use the With function which can be helpful in reducing the size of the code. You con't have to type in all the parameter values. You can use the output from FindMinimum (i.e., sol[[2]]) with ReplaceAll to minimize the typing.

With[
 {
  dataMeas = data[[All, 1 ;; 3]],
  parameterVectors = {c1, b1, c2, b2, cj, bj} /. sol[[2]]
  },
 dataRec = Map[Function[fingers,
    {
     ω1rec[Sequence @@ parameterVectors, Sequence @@ fingers],
     ω2rec[Sequence @@ parameterVectors, Sequence @@ fingers],
       Jrec[Sequence @@ parameterVectors, Sequence @@ fingers]
     }],
   dataMeas
   ]
 ]

Here is a table of the reconstructed (left) and measured data (right).

Now let's plot it.

Column[{
  Show[
   ListPlot[
    Transpose[{
      data[[All, 4]],
      dataRec[[All, 1]]
      }],
    PlotStyle -> {PointSize[Large], Red},
    ImageSize -> 250,
    PlotRange -> {{7, 9}, {7, 9}},
    AxesLabel -> Style["Rec ω1", Bold, 12]
    ],
   Plot[x, {x, 7, 9}, PlotStyle -> {Thin, Black}]
   ],
  Spacer[{1, 10}],
  Show[
   ListPlot[
    Transpose[{
      data[[All, 5]],
      dataRec[[All, 2]]
      }],
    PlotStyle -> {PointSize[Large], Red},
    ImageSize -> 250,
    PlotRange -> {{6.5, 8.5}, {6.5, 8.5}},
    AxesLabel -> Style["Rec ω2", Bold, 12]
    ],
   Plot[x, {x, 6.5, 8.5}, PlotStyle -> {Thin, Black}]
   ],
  Spacer[{1, 10}],
  Show[
   ListPlot[
    Transpose[{
      data[[All, 6]],
      dataRec[[All, 3]]
      }],
    PlotStyle -> {PointSize[Large], Red},
    ImageSize -> 250,
    PlotRange -> {{0.1, 0.6}, {0.1, 0.6}},
    AxesLabel -> Style["Rec J", Bold, 12]
    ],
   Plot[x, {x, 0.1, 0.6}, PlotStyle -> {Thin, Black}]
   ]
  }]

The results look good:

Until new measurements are made we can model the capacitors as:

C1[F1_] := 0.000341*F1 + 0.001572
C2[F2_] := 0.0003019*F2 + 0.003162
Cj[Fj_] := 0.0001930*Fj - 0.00004236

Answer 3

I can unfortunately not just comment because I lack the reputation, which I would have found more appropriate. This case, to me, seems more like a question of actual fitting than the programing so I'll give you a few ideas that I had.

First of all, it could be useful to have analytic expressions for the invcmat elements. This has a simple reason. If those expressions are symmetric in C₁, C₂ your measured data is not giving you enough information to say which of your "identical" condensators has which value, you can just say that one of them has one value while the other one has the other - if the fits converge at all. This would mean that you need an additional measurement to tell them apart and then know which fitted value in each case belongs to which of those.

Second, my fitting routine would be to try a simple least square of differences fit simultaneously fitting \omega₁, \omega₂ and J with the three capacities as parameters (I know that professional software like origin can do this, I am unsure about mathematicas capacities in this point). This could be done for every measurement and then the capacities could be plotted over the number of fingers and you could then think about which function could be appropriate for the look of the data points (separately for each condensator). However, it should be noted that with 4/5 points a lot of functions will yield a graphically pleasing fit, so without any kind of theoretical model, this could end in a bit of guesswork.

I hope I did not misunderstand your question and thus oversimplify it. If I did so, I am sorry in advance for wasting your time.