Thread Subject: simple parameter estimation

please, dear friends,

how'd you solve a system like this:

for a given y1, y2 and y3:

dy1/dt = a*y1 + b*y2 - c*y3
dy2/dt = b*y1 - (a-1)*y2 - c*y1*y3
dy3/dt = a*y1 - b

How to find a, b and c?

Many thanks

Please, is there a way to solve this kind of system by using a ode solver?

Many thanks!

"Marcio " <marciobarbalho@yahoo.com> wrote in message <gcm6me$o7v$1@fred.mathworks.com>...
> Please, is there a way to solve this kind of system by using a ode solver?
>
> Many thanks!

Suppose you knew the parameters in question?
Could you solve the system, predicting y1, y2,
y3? Yes, since you would use an appropriate ode
solver, like ode45, etc.

So what you do is to use a nonlinear regression
tool from the optimization toolbox, to find the
set of parameters that give you the best set of
parameters, a,b,c, that will fit your curves.

Make sure you choose a good set of starting
values for the parameters.

John

> Suppose you knew the parameters in question?
> Could you solve the system, predicting y1, y2,
> y3? Yes, since you would use an appropriate ode
> solver, like ode45, etc.

I don't know the parameters in question. I did not understand what you said, I am trying to find a, b and c, all I know are the "responses", y1, y2 and y3. So, how'd you use ode45 to find the parameters a, b and c?
of course I could give a starting guess if necessary, or a typical solution, as they say!
Many thanks in advance!

"Marcio " <marciobarbalho@yahoo.com> wrote in message <gcmc12$514$1@fred.mathworks.com>...
> > Suppose you knew the parameters in question?
> > Could you solve the system, predicting y1, y2,
> > y3? Yes, since you would use an appropriate ode
> > solver, like ode45, etc.
>
> I don't know the parameters in question. I did not understand what you said, I am trying to find a, b and c, all I know are the "responses", y1, y2 and y3. So, how'd you use ode45 to find the parameters a, b and c?
> of course I could give a starting guess if necessary, or a typical solution, as they say!
> Many thanks in advance!

Then I guess it is time for you to learn what
nonlinear regression is. If you have the
optimization toolbox,

help lsqnonlin
help lsqcurvefit

John

"John D'Errico" <woodchips@rochester.rr.com> wrote
> Then I guess it is time for you to learn what
> nonlinear regression is. If you have the
> optimization toolbox,
>
> help lsqnonlin
> help lsqcurvefit
>
> John

I guess those functions work out with non-linear equations, but it's with the derivatives that I am worried about. How to tell MATLAB they are not just non-linear equations?

Thanks!

> "John D'Errico" <woodchips@rochester.rr.com> wrote
> > Then I guess it is time for you to learn what
> > nonlinear regression is. If you have the
> > optimization toolbox,
> >
> > help lsqnonlin
> > help lsqcurvefit
> >
> > John
>
> I guess those functions work out with non-linear
> equations, but it's with the derivatives that I am
> worried about. How to tell MATLAB they are not just
> non-linear equations?
>
> Thanks!

Assume for a moment that you know the analytical
solution for your system of ordinary differential
equations, say y_i(t,a,b,c) (i=1,2,3).

Then the problem to determine optimal parameters
can be formulated as

min_{a,b,c}:
sum_{j=1}^{n} sum_{i=1}^{3}
(y_i(t_j,a,b,c)-y_measured_i(t_j))^2

where the t_j are times where measurements for the
y_i are available.

This would define a normal nonlinear
parameter fitting problem for which you could use
lsqnonlin, lsqcurvefit.

Now in your case an analytical solution for the
above system of differential equations
seems hard to determine.

For this reason, you have to couple ODE45 and one of the
parameter estimation routines:

In each step k, you get suggestions for parameters
a_k, b_k, c_k from lsqnonlin or lsqcurvefit.
With these parameters, you call ODE45 to
determine the y_i(t_j,a_k,b_k,c_k) and transfer
these values to lsqnonlin or lsqcurvefit.
In the next step, lsqnonlin or lsqcurvefit will
transfer new a_(k+1), b_(k+1), c_(k+1).
You call ODE45 and so on.

I don't know how this can be implemented in MATLAB
in practise, but that's the general idea how you
will have to proceed.

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message
> Now in your case an analytical solution for the
> above system of differential equations
> seems hard to determine.
>
> For this reason, you have to couple ODE45 and one of the
> parameter estimation routines:
>
> In each step k, you get suggestions for parameters
> a_k, b_k, c_k from lsqnonlin or lsqcurvefit.
> With these parameters, you call ODE45 to
> determine the y_i(t_j,a_k,b_k,c_k) and transfer
> these values to lsqnonlin or lsqcurvefit.
> In the next step, lsqnonlin or lsqcurvefit will
> transfer new a_(k+1), b_(k+1), c_(k+1).
> You call ODE45 and so on.

Actually the problem I am trying to solve is pretty much harder than that, I guess it's about over 250 equations, I mean, 250 differential equations and 248 algebraic equations, 598 variables (of course), which gives us a differential-algebraic system of equations, which is, in that case, strictly non-linear.
Of course, I am trying to solve a simple ODE system before starting with a DAE system.

Ok, let me see whether I got your suggestion straight...
You said to solve the system as albebraic and the result feeds the differential system and so on. I think I misunderstood something, will the results from the differential system feed a new algebraic system?

Thanks for commenting!

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in message
> > Now in your case an analytical solution for the
> > above system of differential equations
> > seems hard to determine.
> >
> > For this reason, you have to couple ODE45 and one
> of the
> > parameter estimation routines:
> >
> > In each step k, you get suggestions for parameters
> > a_k, b_k, c_k from lsqnonlin or lsqcurvefit.
> > With these parameters, you call ODE45 to
> > determine the y_i(t_j,a_k,b_k,c_k) and transfer
> > these values to lsqnonlin or lsqcurvefit.
> > In the next step, lsqnonlin or lsqcurvefit will
> > transfer new a_(k+1), b_(k+1), c_(k+1).
> > You call ODE45 and so on.
>
> Actually the problem I am trying to solve is pretty
> much harder than that, I guess it's about over 250
> equations, I mean, 250 differential equations and 248
> algebraic equations, 598 variables (of course), which
> gives us a differential-algebraic system of
> equations, which is, in that case, strictly
> non-linear.
> Of course, I am trying to solve a simple ODE system
> before starting with a DAE system.
>
> Ok, let me see whether I got your suggestion
> straight...
> You said to solve the system as albebraic and the
> result feeds the differential system and so on. I
> think I misunderstood something, will the results
> from the differential system feed a new algebraic
> system?
>
> Thanks for commenting!


You should take a look at

http://num.math.uni-bayreuth.de/~kschittkowski/dyn_sys_book.htm

Here you get a reference to a book that will
hopefully answer your questions.
 
A software-package called EASY-FIT is also
available if the implementation in MATLAB is
too time-consuming.

Best wishes
Torsten.

> You should take a look at
>
> http://num.math.uni-bayreuth.de/~kschittkowski/dyn_sys_book.htm
>
> Here you get a reference to a book that will
> hopefully answer your questions.
>
> A software-package called EASY-FIT is also
> available if the implementation in MATLAB is
> too time-consuming.
>
> Best wishes
> Torsten.

I wish I could, the project I work with works out on MATLAB and Fortran Compiler suite, someone told me I could do that by using some "special" fortran subroutines, it is much faster than MATLAB, but I can't use fortran, that's the problem!!!
what I am really worried about is how to solve such DAE system on MATLAB and then, perform the parameter estimation, another big problem.
Thanks for commenting!

To estimate the parameters in any system, you must have an initial guess for the parameters, x0. Therefore, your overall procedure will be iterative as follows:

1. Guess x0.
2. Solve the system of equations with x0.
3. Compute the cost function (e.g. sum of the square of residuals) by comparing the solutions with your data.
4. Use an optimization algorithm (e.g. fminunc, fminsearch) to automate the reptition of steps 1-3 until the parameter estimation is complete to some predefined tolerance.

"Marcio " <marciobarbalho@yahoo.com> wrote in message <gcl2c0$b9$1@fred.mathworks.com>...
> please, dear friends,
>
> how'd you solve a system like this:
>
> for a given y1, y2 and y3:
>
> dy1/dt = a*y1 + b*y2 - c*y3
> dy2/dt = b*y1 - (a-1)*y2 - c*y1*y3
> dy3/dt = a*y1 - b
>
> How to find a, b and c?
>

Do you know the cauchy data? Or is it a free parameters?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote

> > how'd you solve a system like this:
> >
> > for a given y1, y2 and y3:
> >
> > dy1/dt = a*y1 + b*y2 - c*y3
> > dy2/dt = b*y1 - (a-1)*y2 - c*y1*y3
> > dy3/dt = a*y1 - b
> >
> > How to find a, b and c?
> >
>
> Do you know the cauchy data? Or is it a free parameters?

t = [t0 t1 t2];
y1 = [x1 x2 x3];
y2 = [p1 p2 p3];
y3 = [w1 w2 w3];

I've been trying to transform that system into a simple non-linear equation system by transforming the derivatives into new functions, for example, dy1/dt could be called f1(t) once y1 is known. I'm working on it. I just had this insight, whether it works out, my day will be made!

"Marcio " <marciobarbalho@yahoo.com> wrote in message <gcog9f$219$1@fred.mathworks.com>...

> >
> > Do you know the cauchy data? Or is it a free parameters?
>
> t = [t0 t1 t2];
> y1 = [x1 x2 x3];
> y2 = [p1 p2 p3];
> y3 = [w1 w2 w3];
>

OK, so you have state parameters (y) at three times. So your problem is estimation parameters from state data. This problem is well know.

I'll repeat myself again: See papers from Jacques Louis Lions on Control Optimal theory.

What you need is compute the derivative of y (in particular at {t0,t1,t2}) with respect to {a,b,c}. This can be accomplished in two steps:
1. Solving the linearization of the direct ODE
2. Solving an adjoint equation of the linearized ODE, with a forcing term as the deviation to data.

More details on the above reference.

After the derivative is computed, plug it in an optimizer and let it does the works.

Bruno

> I've been trying to transform that system into a simple non-linear equation system by transforming the derivatives into new functions, for example, dy1/dt could be called f1(t) once y1 is known. I'm working on it. I just had this insight, whether it works out, my day will be made!

On Oct 9, 9:53 am, "Marcio " <marciobarba...@yahoo.com> wrote:
> please, dear friends,
>
> how'd you solve a system like this:
>
> for a given y1, y2 and y3:
>
> dy1/dt = a*y1 + b*y2 - c*y3
> dy2/dt = b*y1 - (a-1)*y2 - c*y1*y3
> dy3/dt = a*y1 - b
>
> How to find a, b and c?
>
> Many thanks

If y1,y2,y3 are given, dy1/dt,dy2/dt and dy3/dt
are readily estimated.

You then have 3 systems of N linear equations for
A = [a b c]'.

Define

Y1 = [y1 y2 -y3];
Y2 = [-y2 y1 -y1.*y3];
Y3 = [y1 -ones(N,1) zeros(N,1)];

Then

dy1/dt = Y1*A
dy2 - y2 = Y2*A
dy3/dt = Y3*A

Note that Y3 is rank deficient. Therefore, the LS
solution of the 3rd system may not be unique.
Although Y1 and Y2 are full rank, the 1st two equations
yield 2 unique, but very different LS solutions
for A. In general, all of the three LS solutions will
be different.

A more satisfactory approach is to use a combined LS
solution. If there is evidence that points or equations
have unequal importance, weighting is easily introduced.

The ordinary LS solution results from minimizing


f = ||dy1/dt-Y1*A||^2 + ||dy2/dt-y2-Y2*A||^2
     + ||dy3/dt-Y3*A||^2.


Hope this helps.

Greg

On Oct 10, 11:25 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Oct 9, 9:53 am, "Marcio " <marciobarba...@yahoo.com> wrote:
>
> > please, dear friends,
>
> > how'd you solve a system like this:
>
> > for a given y1, y2 and y3:
>
> > dy1/dt = a*y1 + b*y2 - c*y3
> > dy2/dt = b*y1 - (a-1)*y2 - c*y1*y3
> > dy3/dt = a*y1 - b
>
> > How to find a, b and c?
>
> > Many thanks
>
> If y1,y2,y3 are given, dy1/dt,dy2/dt and dy3/dt
> are readily estimated.
>
> You then have 3 systems of N linear equations for
> A = [a b c]'.
>
> Define
>
> Y1 = [y1 y2 -y3];
> Y2 = [-y2 y1 -y1.*y3];
> Y3 = [y1 -ones(N,1) zeros(N,1)];
>
> Then
>
> dy1/dt = Y1*A
> dy2 - y2 = Y2*A
> dy3/dt = Y3*A
>
> Note that Y3 is rank deficient. Therefore, the LS
> solution of the 3rd system may not be unique.
> Although Y1 and Y2 are full rank, the 1st two equations
> yield 2 unique, but very different LS solutions
> for A. In general, all of the three LS solutions will
> be different.
>
> A more satisfactory approach is to use a combined LS
> solution. If there is evidence that points or equations
> have unequal importance, weighting is easily introduced.
>
> The ordinary LS solution results from minimizing
>
> f = ||dy1/dt-Y1*A||^2 + ||dy2/dt-y2-Y2*A||^2
> + ||dy3/dt-Y3*A||^2.

I forgot how to take the gradient of f w.r.t. a,b,c.
However, it looks like the result should be the same as
if you just added the three equations together.
Therefore, a weighted version would look like

w1*dy1/dt + w2*(dy2-y2) + w3*dy3/dt

   = (w1*Y1 + w2*Y2 + w3*Y3)*A

with solution

A = (w1*Y1 + w2*Y2 + w3*Y3) ...
     \(w1*dy1/dt + w2*(dy2-y2) + w3*dy3/dt)

Hope this helps.

Greg

On Oct 12, 12:36 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Oct 10, 11:25 pm,Greg Heath<he...@alumni.brown.edu> wrote:
>
> > On Oct 9, 9:53 am, "Marcio " <marciobarba...@yahoo.com> wrote:
>
> > > please, dear friends,
>
> > > how'd you solve a system like this:
>
> > > for a given y1, y2 and y3:
>
> > > dy1/dt = a*y1 + b*y2 - c*y3
> > > dy2/dt = b*y1 - (a-1)*y2 - c*y1*y3
> > > dy3/dt = a*y1 - b
>
> > > How to find a, b and c?
>
> > > Many thanks
>
> > If y1,y2,y3 are given, dy1/dt,dy2/dt and dy3/dt
> > are readily estimated.
>
> > You then have 3 systems of N linear equations for
> > A = [a b c]'.
>
> > Define
>
> > Y1 = [y1 y2 -y3];
> > Y2 = [-y2 y1 -y1.*y3];
> > Y3 = [y1 -ones(N,1) zeros(N,1)];
>
> > Then
>
> > dy1/dt = Y1*A
> > dy2 - y2 = Y2*A
> > dy3/dt = Y3*A
>
> > Note that Y3 is rank deficient. Therefore, the LS
> > solution of the 3rd system may not be unique.
> > Although Y1 and Y2 are full rank, the 1st two equations
> > yield 2 unique, but very different LS solutions
> > for A. In general, all of the three LS solutions will
> > be different.
>
> > A more satisfactory approach is to use a combined LS
> > solution. If there is evidence that points or equations
> > have unequal importance, weighting is easily introduced.
>
> > The ordinary LS solution results from minimizing
>
> > f = ||dy1/dt-Y1*A||^2 + ||dy2/dt-y2-Y2*A||^2
> > + ||dy3/dt-Y3*A||^2.
>
> I forgot how to take the gradient of f w.r.t. a,b,c.
> However, it looks like the result should be the same as
> if you just added the three equations together.
> Therefore, a weighted version would look like
>
> w1*dy1/dt + w2*(dy2-y2) + w3*dy3/dt
>
> = (w1*Y1 + w2*Y2 + w3*Y3)*A
>
> with solution
>
> A = (w1*Y1 + w2*Y2 + w3*Y3) ...
> \(w1*dy1/dt + w2*(dy2-y2) + w3*dy3/dt)

My mathematician brother helped me with the
gradient. Tthe resulting equation is

 (w1*Y1^T*dy1/dt + w2*Y2^T*(dy2-y2) + w3*Y3^T*dy3/dt

    = (w1*Y1^T*Y1 + w2*Y2^T*Y2 + w3*Y3^T*Y3)*A

which is readily solved given two ratios of (w1,w2,w3).

However, since w = [w1,w2,w3] is arbitrary, this isn't
the full story.

If w1=w2=w3, then the three individual squared
residuals will probably have different sizes. Therefore,
the three equations really won't have equal importance.
Consequently, one could try to choose w to equalize the
squared residuals!

Since this probably falls into the category of Quadratic
Optimization with quadratic constraints, I'd reformulate
the problem into one of Quadratic Optimization with
Linear constraints.

I guess the following iterative linear approach is also
a possibility:

1. Initialize the permanent value w1 = 1,
2. Initialize w2(0) = 1, w3(0))=1
2. Solve for A(0)
3. Calculate MSE1(0) = ||Y1*A(0) -dy1/dt||^2, etc
4. Update w2(1) = MSE1(0)/MSE2(0), etc.
5. Solve for A(1), etc

Hope this helps.

Greg-

This is a simple parameter estimation problem
you must solve two key tasks:
1. ODE Solving (Using ODE solver or any tool) with given parameters (a =1, b=2, c=3 for example)
2. Optimal parameters searching (Using Optimization tool to find out optimal results, almost tools based on GA, SA.. to find out global optimal results)
So the input of your problem must be observed data of yk(t) (k=1,2,3), t = t1, t2, t3... tn (n observed samples )

You start with initial a, b,c. Then use ODE to find out y1, y2, y3 from equations, compare with observed data to get the fitness value (least square error for example) for your optimization search

After all, you can use many tools (Solver in Excel from Front line.., Matlab tool or just write a Fotran program by yourself with large problem (Visual compact Fortran also has ODE...)). Details just read Matlab help! Nice.