This is similar to this question where the OP had a coupled system of 4 ODEs and subject to 4 data sets could fit the parameters.
While the question wasn't solved directly, Oleksandr R. gave a very good and detailed response on how it could be done.
I have a similar case where I have a system of 4 coupled PDEs. Provided the values I can solve the equations no problem, however given some data and working backwards I am slightly stuck, probably due to the nature of the fitting.
I have 4 PDEs with functions y1[x,t], y2[x,t], z1[x,t], z2[x,t]
ClearAll[y1, y2, z1, z2, a, b, \[Alpha]1, \[Alpha]2, K1, K2, K31, \
K32, y1in, y2in, L, x, t]
y1in = 1;
y2in = 1;
L = 20;
System = {a D[y1[x, t], t] + D[y1[x, t], x] == b D[y1[x, t], x, x] - D[z1[x, t], t],
a D[y2[x, t], t] + D[y2[x, t], x] == b D[y2[x, t], x, x] - \[Alpha]2/\[Alpha]1 y1in/y2in D[z2[x, t],t],
D[z1[x, t], t] == y1[x, t] (1 - z1[x, t] - z2[x, t]) - 1/(K1 y1in) z1[x, t] - y2in/y1in * 1/K31 y2[x, t] z1[x, t],
y2in/y1in * 1/K31 D[z2[x, t], t] == y2[x, t] (1 - z1[x, t] - z2[x, t]) - 1/(K2 y2in) z2[x, t] + K31/K32 * 1/K31 y2[x, t] z1[x, t],
z1[x, 0] == 0,
z2[x, 0] == 0,
(y1[x, 0]) == If[x == 0, y1in, 0],
(y2[x, 0]) == If[x == 0, y2in, 0],
y1[0, t] == y1in,
y2[0, t] == y2in,
D[y1[x, t], x] == 0 /. x -> L,
D[y2[x, t], x] == 0 /. x -> L};
When solving with known parameter values I had a lot of success using MethodOfLines in NDSolveValue. The solutions match experimental data well so I decided to keep the method for the ParametricNDSolveValue
soln = ParametricNDSolveValue[
System, {y1, y2, z1, z2}, {x, 0, L}, {t, 0, 250}, {a,
b, \[Alpha]1, \[Alpha]2, K1, K2, K31, K32},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 101}}];
It runs with no errors and I get interpolating functions of parameters. Here is where my question differs from linked question.
I have data for y1 and y2 only so we are restricted to only fitting y1 and y2 while still needing to solve for unknown parameters in z1 and z2. The experimental data I have was taken at x=L so my data is also restricted to being evaluated at a given x. It is a set of time and function values at x=L and at those times. {ti,y1[x=L,ti]}.
y1data = {{0, 0}, {5.53333332762122, 0}, {11.9166666653473,
0}, {18.3166666631587, 0}, {24.7166666609701,
0}, {31.1166666587815, 0}, {37.5166666670702,
0}, {43.8999999943189, 0}, {50.283333332045,
0}, {56.6666666592937, 0}, {63.0666666675825,
0}, {69.4333333347458, 0}, {75.8166666619945,
0}, {82.1333333279472, 0}, {88.5333333257586,
0}, {94.9333333340474, 0}, {101.333333331859,
0}, {107.733333329670, 0}, {114.133333327482,
0}, {120.533333325293, 0}, {126.933333333582,
0}, {133.333333331393, 0}, {139.716666658642,
0}, {146.099999996368, 0}, {152.483333334094,
0.0254801793003281}, {158.883333331905,
0.0570225961819381}, {165.266666659154,
0.139066269334469}, {171.64999999688,
0.275847461164117}, {178.049999994691,
0.455659062262623}, {184.449999992503,
0.603469144149784}, {190.850000000792,
0.746343875992287}, {197.23333332804,
0.913864055626508}, {203.599999995204,
0.974246446130602}, {209.999999993015,
0.994394777919682}, {216.33333332953,
0.999830966869506}, {222.666666666046,
0.997514586933113}, {228.966666661436,
0.999361430395913}, {235.283333327388,
0.999997913171228}, {241.666666665114, 0.999528376697635}};
y2data = {{0, 0}, {5.53333332762122, 0}, {11.9166666653473,
0}, {18.3166666631587, 0}, {24.7166666609701,
0}, {31.1166666587815, 0}, {37.5166666670702,
0}, {43.8999999943189, 0}, {50.283333332045,
0.0002}, {56.6666666592937,
0.000230340266319416}, {63.0666666675825,
0.48369536424858}, {69.4333333347458,
1.45053583544287}, {75.8166666619945,
1.47667305732939}, {82.1333333279472,
1.48662247716624}, {88.5333333257586,
1.48668646057355}, {94.9333333340474,
1.50270790576421}, {101.333333331859,
1.47977625258397}, {107.73333332967,
1.48325694994169}, {114.133333327482,
1.49327035318585}, {120.533333325293,
1.49996941593131}, {126.933333333582,
1.47763280843905}, {133.333333331393,
1.49058305007879}, {139.716666658642,
1.49380141546653}, {146.099999996368,
1.48882350637774}, {152.483333334094,
1.48003218621321}, {158.883333331905,
1.44971044948855}, {165.266666659154,
1.43846216648329}, {171.64999999688,
1.32757892161342}, {178.049999994691,
1.26303885865886}, {184.449999992503,
1.18215103513636}, {190.850000000792,
1.08591999054069}, {197.23333332804,
1.01503277358089}, {203.599999995204,
1.01980593576629}, {209.999999993015,
1.01160966128976}, {216.33333332953,
0.991704423275321}, {222.666666666046,
1.01802719704305}, {228.966666661436,
0.989631360878446}, {235.283333327388,
1.00947261548557}, {241.666666665114, 0.997347759800146}};
Can I still solve for all my parameters given these 2 data sets? The PDEs are linked so theoretically it should be possible yes?
I am assuming I need to go through a similar procedure of setting up the abscissae but now for time only as I need to evaluate the interpolating functions at x=L. But how do I specify that the solutions for y1 and y2 must fit the data for y1 and y2 while still calculating the value of the parameters from the given data?