 |
 | Home | Contact Us | FAQ | Search & Site Map | Link to Us |
 | Sign In | Join | Other 45 Sites in Network |
Math Forum / Math Software / MATLAB / July 2008Solving an inverse ODE problem
Hi, I have a problem dealing with the following ODE system: C * dx(t)/dt = A - B * x(t) - y(t) where A,B and C are constants. For this system I have data from measurements giving a noisy, continous signal of x. My goal is to get an continous signal for y(t). By the use of splines a first estimate of dx/dt(t) and x(t) can be achieved at any time resulting in an estimate of y(t). To check whether y(t) is calculated correctly the upper equation is solved using the estimated y(t). Comparing the results of the calculation with the measurement gives significant deviations at some points which indicates a wrong estimate of y(t). Is there a method to adapt y(t) in a certain way, that the next calculation will lead to a better agreement of calculation and measurement? The idea is to run this optimization as long as the measurements and the calculated data agree adequately. This would mean that y(t) is calculated correctly. The main problem is, that C is quite big, so any changes in y(t) and dx(t)/dt will result in an corresponding change of x with a certain delay. > Hi, > [quoted text clipped - 36 lines] > change of > x with a certain delay. Let t_0 < t_1 < t_2 <...< t_n be times where measurements x_measured(t_i) are available. Let y_i~ be the value for y~=y/C on the intervall [t_(i-1),t_i]. Then solve the following problem: min: sum_{i=1}^{n} (x(t_i)-x_measured(t_i))^2 where x is given by the differential equation dx/dt = A~ - B~*x(t) - y_i~ if t_(i-1) <= t <= t_i. So x(t_i) depends on the parameters A~, B~, y_1~ ,..., y_i~ which have to be calculated by nonlinear regression. (The whole set of unknown parameters which have to be calculated are A~, B~, y_1~,...,y_n~). Best wishes Torsten. > > Hi, > > [quoted text clipped - 72 lines] > Best wishes > Torsten. To make the solution y~ "more continuous", add a parameter y_0~ at time t0. Then, given iterates y_0~,y_1~,...,y_n~ for the unknown values y_i~, you could first approximate them by a spline y~_spline. Then the problem to solve looks like min: sum_{i=1}^{n} (x(t_i)-x_measured(t_i))^2 where x is given by the differential equation dx/dt = A~ - B~*x(t) - y~_spline(t). In this way, the discontinuities in the x-derivative at times t_i have disappeared. Best wishes Torsten. Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message > So x(t_i) depends on the parameters A~, B~, y_1~ ,..., > y_i~ which have to be calculated by nonlinear regression. I would like to add a comment. As A, B, C are constant and known in the problem, x(t_i) depends even linearly to y. The forcing y can be thus estimated by *linear* regression, with all the nice properties coming in bonus: minimum exists, no local minima, etc... (the minimum is even unique IMO in this example). Several reference text books of this type of problem is from Jacques Louis Lions http://www.ann.jussieu.fr/jllions/. I'm a great admirer of him. Weather forecast and climatology modeling own a big debt to him. Bruno "Bruno Luong" <b.luong@fogale.fr> wrote in message <g6apii$pfs$1@fred.mathworks.com>... > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message > > [quoted text clipped - 15 lines] > > Bruno Thanks for your answers, A, B and C are indeed known constants in my case. Can you recommend specific reference? I couldn't really figure out what is suitable for my problem. S H "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message <g6brlr$50h$1@fred.mathworks.com>... > A, B and C are indeed known constants in my case. Can you > recommend specific reference? I couldn't really figure out > what is suitable for my problem. > > S H 1. Jacquess Louis Lions: CONTROLE OPTIMAL DE SYSTEMES GOUVERNES PAR DES EQUATIONS AUX DERIVEES (Book in french) This is a reference text book. That might not explain the practical aspects, but gives a solid background and framework on Adjoint equation and Optimality System. 2. F.-X. Le Dimetet O. Talagrand. Variational algorithms for analysis and assimilation of meteorogical observations : theoretical aspects. Tellus, 38A: 97-110,1986. This is a basic paper for assimilation. The equation is PDE with a general form, but many steps are described in details. Bruno > "Bruno Luong" <b.luong@fogale.fr> wrote in message > <g6apii$pfs$1@fred.mathworks.com>... [quoted text clipped - 40 lines] > > S H The problem can be solved as follows: You have discrete times t_0 < t_1 < t_2 <...< t_n where you have measurements for your function x: x_measured(t_0),...,x_measured(t_n). If t_i and t_(i+1) are "close", we may assume that the "control variable" y is constant on each of the intervals [t_i,t_(i+1)] and equals y_(i+1). Then, solving your differential equation over the time intervall [t_i,t_i+1] gives log ((A/C+B/C*x(t_(i+1))+y_(i+1)/C)/ (A/C+B/C*x(t_i)+y_(i+1)/C)) = t_(i+1)-t_i or equivalently x(t_(i+1))=(A/B+x(t_i)+y_(i+1)/B))*e^(t_(i+1)-t_i) -A/B-y_(i+1)/B (i=0,...,n-1) This is the simulated solution for x given the control vector [y_1,...,y_n]. Now you can simply use a linear regression tool to solve the following regression problem in the unknown parameter vector [y_1,...,y_n]: minimize: sum_{i=1}^{n} (x(t_(i))-x_measured(t_i))^2 Of course, if you have constraints on the vector [y_1,...,y_n] (e.g. y(t) <=5 or y(t) >= -0.5), the problem will become more difficult. Best wishes Torsten. > The problem can be solved as follows: > [quoted text clipped - 16 lines] > > (i=0,...,n-1) Sorry, I forgot the B/C as inner derivative: .. log ((A/C+B/C*x(t_(i+1))+y_(i+1)/C)/ (A/C+B/C*x(t_i)+y_(i+1)/C)) = B/C*(t_(i+1)-t_i) or equivalently x(t_(i+1))=(A/B+x(t_i)+y_(i+1)/B))*e^(B/C*(t_(i+1)-t_i)) -A/B-y_(i+1)/B (i=0,...,n-1) .. > This is the simulated solution for x given the > control vector [y_1,...,y_n]. [quoted text clipped - 11 lines] > Best wishes > Torsten. Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <1389485.1216972746386.JavaMail.jakarta@nitrogen.mathforum.org>... > > The problem can be solved as follows: > > [quoted text clipped - 45 lines] > > Best wishes > > Torsten. Thanks for the useful hints, there's another question following from the solution of the differential equation. Integrations rules tell me that the ODE is solved as log |((A/C-B/C*x(t_(i+1))-y_(i+1)/C)/ (A/C-B/C*x(t_i)-y_(i+1)/C))| = B/C*(t_(i+1)-t_i) with using some "-" instead of "+" and using the absolute value (| ... |). Is there any reason you derived it differently? Is there a MATLAB solver which can handle my problem? S H Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <1389485.1216972746386.JavaMail.jakarta@nitrogen.mathforum.org>... > > The problem can be solved as follows: > > [quoted text clipped - 45 lines] > > Best wishes > > Torsten. Thanks for the useful hints, there's another question following from the solution of the differential equation. Integration rules tell me that the ODE is solved as log |((A/C-B/C*x(t_(i+1))-y_(i+1)/C)/ (A/C-B/C*x(t_i)-y_(i+1)/C))| = B/C*(t_(i+1)-t_i) with using some "-" instead of "+" and using the absolute value (| ... |). Is there any reason you derived it differently? Is there a MATLAB solver which can handle my problem? S H >Thanks for the useful hints, > [quoted text clipped - 8 lines] >value (| ... |). Is there any reason you derived it >differently? No, I just had in mind that your ODE reads C*dx/dt = A + B*x + y So the solution is given by log(A-B*x(t_(i+1))-y(t_(i+1)))/(A-B*x(t_i)-y(t_(i+1)))) = -B/C*(t_(i+1)-t_i) or equivalently x(t_(i+1)) = A/B - y(t_(i+1))/B - 1/B*(A-B*x(t_i)-y(t_i+1)))*e^(-B/C*(t_(i+1)-t_i)) (i=0,...,n-1) So you can recursively compute your unknowns y(t_(i+1)), x(t_(i+1)) as follows: For i = 0,...,n-1 do Calculate y(t_(i+1)) from A/B - y(t_(i+1))/B - 1/B*(A-B*x(t_i)-y(t_(i+1)))*e^(-B/C*(t_(i+1)-t_i)) - x_measured(t_(i+1)) = 0 (linear equation in y(t_(i+1))) Calculate x(t_(i+1)) from x(t_(i+1)) = A/B - y(t_(i+1))/B - 1/B*(A-B*x(t_i)-y(t_i+1)))*e^(-B/C*(t_(i+1)-t_i)) end for >Is there a MATLAB solver which can handle my >problem? > >S H Best wishes Torsten. Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <13474757.1217241476173.JavaMail.jakarta@nitrogen.mathforum.org>... > >Thanks for the useful hints, > > [quoted text clipped - 44 lines] > Best wishes > Torsten. Ok, do you recommend a certain tool/solver? S H >Ok, > >do you recommend a certain tool/solver? > >S H In FORTRAN, it is one loop: Save the times t_i for which you have measurements in an array t(i) Save the measured x at time t_i in x_measured(i) Give values to A, B and C x(0) = xstart Do i = 0, n-1 y(i+1) = (A-(A-B*x(i))*exp(-B/C*(t(i+1)-t(i)))- x_measured(i+1)*B)/(1-exp(-B/C*(t(i+1)-t(i))) x(i+1) = A/B - y(i+1)/B - 1/B*(A-B*x(i)-y(i+1))* exp(-B/C*(t(i+1)-t(i))) end do The array y(i) (i=1,...,n) contains the control vector at times t(i) I do not work with MATLAB ; so I don't know how such a Do-loop is programmed in this environment. If the measured x_measured(i) contain much noise, it may be advisable first to smoothe these data (e.g. by spline approximation). Best wishes Torsten. Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <33249462.1217257035765.JavaMail.jakarta@nitrogen.mathforum.org>... > If the measured x_measured(i) contain much noise, > it may be advisable first to smoothe these data > (e.g. by spline approximation). A better way is to leave the noise in the data x_measured and perform regularization on the control y by adding a Tchikonov regularization term in the control. In other word instead of solving original least-square problem minimize |A*y - x_measured|^2 Solve the penalized least-square function: minimize |A*y - x_measured|^2 + lambda*|y|^2 lambda is some positive number, that is larger if x_measured is noisier. Bruno "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <g6koaj$sqf$1@fred.mathworks.com>... > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in > message <33249462.1217257035765.JavaMail.jakarta@nitrogen.mathforum.org>... > > If the measured x_measured(i) contain much noise, > > it may be advisable first to smoothe these data [quoted text clipped - 15 lines] > > Bruno Ok, but back to the integrated function. The integration of "dx/(ax+c)" is "1/a log |ax+c|". I don't see why Torsten simply neglected the bars for the absolute values. In minimize |A*y - x_measured|^2: what is the expession for A? S H "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message <g6n3jt$69r$1@fred.mathworks.com>... > In minimize |A*y - x_measured|^2: what is the expession for A? A represents the matrix that maps linearly the control "y" to your solution x taken at measuring points. Bruno "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <g6n4nh$i0b$1@fred.mathworks.com>... > "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message > <g6n3jt$69r$1@fred.mathworks.com>... [quoted text clipped - 5 lines] > > Bruno Yea, but still my problem with the absolute values exist, so I don't see how to get the A. "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message <g6nah5$t07$1@fred.mathworks.com>... > Yea, but still my problem with the absolute values exist, so > I don't see how to get the A. When you solve formally the ODE dx/dt = ax + b (Nore: consider a, b as function of your parameter A, B, C, y, so the above eqt is general as yours) You get log | ax + b | = a*t + cte1 OR | ax + b | = cte2 * exp(a*t) But this is the same as ax + b = cte3 * exp(a*t) (just selecting cte3 = +/- cte2 accordingly), OR x = c * exp(a*t) - b/a If you have a cauchy data x0, i.e., x(t=0) = x0, then c must be equal to x0 + b/a So the solution of Cauchy problem is x(t) = (x0 + b/a) * exp(a*t) - b/a With or without absolute sign does not really matter, it's essentially the same solution. You notice that that a*x(t) + b = (x0 + b/a) * exp(a*t) never change the sign for t>0. Bruno I made a mistake in one of the equations, but those who still follow will figure it out anyway. Bruno "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <g6nvqa$5i4$1@fred.mathworks.com>... > I made a mistake in one of the equations, but those who > still follow will figure it out anyway. > > Bruno That's what I thought as well but solving my equation brings me to |(A-B*x(t_i+1)-y) / (A-B*x(t_i)-y) | = exp(-B/C(t_i+1-t_i)) with the right hand side being always postive. When the absolute values are neglected on the left hand side it can get negative e.g. when "B*x(t_i+1) > A-y" and "B*x(t_i+1) < A-y" which may occur when dx/dt changes sign. S H "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message <g6p3ss$8q1$1@fred.mathworks.com>... > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in > message <g6nvqa$5i4$1@fred.mathworks.com>... [quoted text clipped - 14 lines] > > S H Sorry, the critical point is e.g. when "B*x(t_i+1) > A-y" and "B*x(t_i) < A-y" You haven't listen: it must *not* change sign if A and B and y are constant in the interval. If y changes, then probably yes. But then you only have to split the your interval at the place where y suppose to change. That's why Tosten assumes y is constant on each sub-interval. You are free to select to the space of your control and the sampling size and place. You ODE is linear and Lizschizien and all, there is a unique solution for Cauchy problem, there is not two or more separate solutions. The exponential expression provides an analytic solution for Cauchy problem (if you don't believe take the derivative of it and replace to the equation, verity the initial condition). So THE solution is given there on each sub-interval. Bruno > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in > message <g6nvqa$5i4$1@fred.mathworks.com>... [quoted text clipped - 10 lines] > |(A-B*x(t_i+1)-y) / (A-B*x(t_i)-y) | = > exp(-B/C(t_i+1-t_i)) The fact that the right hand side is positive shows that (A-B*x(t_i+1)-y) and (A-B*x(t_i)-y) are both positive or both negative - so absolute values are no longer necessary. > the > absolute values are neglected on the left hand side [quoted text clipped - 4 lines] > > S H Best wishes Torsten. > in > > message <g6nvqa$5i4$1@fred.mathworks.com>... [quoted text clipped - 17 lines] > positive or both negative - so absolute values > are no longer necessary. This argumentation is of course rubbish. But see Bruno's last reply. > > the > > absolute values are neglected on the left hand [quoted text clipped - 8 lines] > Best wishes > Torsten. Best wishes Torsten. > >Ok, > > [quoted text clipped - 32 lines] > Best wishes > Torsten. It's even easier than that: by the integration of the differential equation between t_i and t_(i+1), we got y_(i+1)-y_(i+1)*e^(-B/C*(t_(i+1)-t_i)) = A - B*x_(i+1) - (A-B*x_i)*e^(-B/C*(t_(i+1)-t_i)) where y_(i+1) is the (constant) control on the interval [t_i,t_(i+1)]. By starting with x(t_0) = x_measured(t_0), we can directly choose y_(i+1) such that x_(i+1) equals x_measured(t_(i+1)) by solving y_(i+1)-y_(i+1)*e^(-B/C*(t_(i+1)-t_i)) = A - B*x_measured(t_(i+1)) - (A-B*x_measured(t_i))*e^(-B/C*(t_(i+1)-t_i)) for y_(i+1): y_(i+1) = (A - B*x_measured(t_(i+1)) - (A-B*x_measured(t_i))*e^(-B/C*(t_(i+1)-t_i)))/ (1-e^(-B/C*(t_(i+1)-t_i))) Best wishes Torsten. > > >Ok, > > > [quoted text clipped - 59 lines] > Best wishes > Torsten. I thought about this "easy" solution again and now see one problem: The state x_measured(i+1) may not be reachable by a control y_(i+1) from the state x_measured(i). So you have to stick to the old formulation as a nonlinear regression problem: Minimize : F(y_1,...,y_n) = sum_{i=0}^{n-1} (x_(i+1)-x_measured(t_(i+1)))^2 where the x_(i+1) are given by the solution of your differential equation: x_(i+1) = 1/B*(A-y_(i+1)-(A-B*x_i-y_(i+1))* e^(-B/C*(t_(i+1)-t_i))) for i = 0,...,n-1. starting with x_0 = x_measured(t_0) In my opinion, the MATLAB solver LSQCURVEFIT is the correct choice for this kind of problem. Best wishes Torsten. Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message <24187566.1217406182906.JavaMail.jakarta@nitrogen.mathforum.org>... > > > >Ok, > > > > [quoted text clipped - 78 lines] > e^(-B/C*(t_(i+1)-t_i))) > for i = 0,...,n-1. This recursive formula on {x} is still linear (well affine to be 100% correct) to control {y}. There is no need for non-linear fit. Another remark is the ode solution inherently decreases as exponential with constant time scale. Matching exactly xi to xi_data is possible but that might requires huge variation of rhs yi on noisier data at high sample. Adding penalization on y will implies a smother convergence, providing only approximation of the data. This is very similar to filtering, PID. Bruno > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote > in [quoted text clipped - 108 lines] > > Bruno Two times agreed. So without the penalty term in the objective function of the regression problem, it should be possible to calculate the y_i directly. With the penalty term: which MATLAB solver can be used ? Best wishes Torsten. Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message > With the penalty term: which MATLAB solver can be > used ? I'm not sure MATLAB has a canned solver. There are two approaches (one can show they are sometime equivalent) 1. Statistic approach. For convenient, the ODE after discretiation is x_i = 1/B*(A-y_(i)-(A-B*x_(i-1)-y_(i))* e^(-B/C*(t_(i)-t_(i-1)))). In other word, synthetize it x_i = f_i*x_(i-1) + g_i*y_i + h_i (eq1) with f_i, g_i, h_i are functions of A, B, C, t(i) and t(i-1). If x_i_data = x_i + noise_i, with noise_i is unknown Gaussian white noise with known standard deviation. One can estimate the control y_i using Kalman filter. 2. Deterministic approach. Write (eq1) for i=1,2,... x_1 = f_1 x_0 + g_1 y_0 + h_1 x_2 = f_2*f_1 x_0 + (f_2*g_1)*y_0 + g_2*y_1 + (f_2*h_1 + h_2) ... x_n = ... This can be written in Matrix form x = A*y + b, with x = (x_1,...,x_n)' y = (y_1,...,y_n)' A = some triangular matrix that can be computed recursively by column, and b is some fix vector in R^n, also computes recusively. Least square estimation of y from x_data is done by y = -pinv(A) * b (or = A\b for underdetermined system) = -inv(A'*A) * A' *b Tchikonov penalization is doing this: y = -inv[A'*A + lambda I] * A' *b. Bruno "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <g6qb52$d5i$1@fred.mathworks.com>... > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message > > With the penalty term: which MATLAB solver can be [quoted text clipped - 52 lines] > > Bruno Back to the main problem. I still have difficulties with the recursive character of my equation. I can't really find a suitable solver for it. "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message <g6rqt6$ncp$1@fred.mathworks.com>... > Back to the main problem. I still have difficulties with the > recursive character of my equation. I can't really find a > suitable solver for it. Why do you need a solver? The recursive formula is given explicitly by Torsten, can't you just put it in a matlab code? Bruno "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <g6rudm$e78$1@fred.mathworks.com>... > "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message > <g6rqt6$ncp$1@fred.mathworks.com>... [quoted text clipped - 7 lines] > > Bruno yea, but I need a minimizer, right? > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in > message <g6rudm$e78$1@fred.mathworks.com>... [quoted text clipped - 16 lines] > > yea, but I need a minimizer, right? Explicit solution for the controls is y_(i+1) = (A - B*x_measured(t_(i+1)) - (A - B*x_measured(t_i))* e^(-B/C*(t_(i+1)-t_i)))/ (1-e^(-B/C(t_(i+1)-t_i))) ; no minimizer needed. Best wishes Torsten. "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message <g6rvoq$q32$1@fred.mathworks.com>... > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in > message <g6rudm$e78$1@fred.mathworks.com>... > > "S H" <sven.hansen@avt.rwth-aachen.de> wrote in message > > <g6rqt6$ncp$1@fred.mathworks.com>... > yea, but I need a minimizer, right? Not necessary, if you follows the deterministic approach MATLAB linear algebra functions (such as "\") is OK, see above. If you prefer brute force solution, fmincon would do (why the hell do we need a non-linear tool to solve a linear problem). Matlab possibly has a linear least square routine where user is not required to provide explicitly the matrix. I don't know optimization toolbox to advise you. So far I have always managed to avoid it. Bruno "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <g6s106$7li$1@fred.mathworks.com>... > Matlab possibly has a linear least square routine where user > is not required to provide explicitly the matrix. I don't > know optimization toolbox to advise you. So far I have > always managed to avoid it. I found it, it calls "lsqr". To use it with a fun form you need to program not only the direct formula but also its adjoint. That may be tricky and requires a minimum of skill. Bruno >Back to the main problem. I still have difficulties with >the >recursive character of my equation. I can't really find a >suitable solver for it. Intgrating your equation C*dx/dt = A - B*x - y between t_i and t_(i+1) and assuming constant control y_(i+1) on that intervall gave A-B*x_(i+1)-y_(i+1) = e^(-B/C*(t_(i+1)-t_i))* (A-B*x_i-y_(i+1)) Now you can find controls y_i such that your solution x_i of the differential equation _exactly coincides_ with your measurements in the t_i's by solving the above equation for y_(i+1) and replacing x_i, x_(i+1) by x_measured(t_i), x_measured(t_(i+1)): y_(i+1) = (A - B*x_measured(t_(i+1)) - (A - B*x_measured(t_i))* e^(-B/C*(t_(i+1)-t_i)))/ (1-e^(-B/C(t_(i+1)-t_i))) This is a direct equation to solve for the controls - no recursion necessary. Maybe you should check the results y_i obtained by this method first. If the jumps in the controls from interval to interval are too large (maybe due to the noise in your measurements), it's still time to concentrate on more sophisticated methods suggested by Bruno. Best wishes Torsten. "Bruno Luong" <b.luong@fogale.fr> wrote in message <g6apii$pfs$1@fred.mathworks.com>... > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in message > > [quoted text clipped - 15 lines] > > Bruno Thanks for your answers, A, B and C are indeed known constants in my case. Can you recommend specific reference? I couldn't really figure out what is suitable for my problem. S H |
Reply to this Message |
 Sign In
 Join
 My Latest Posts My Monitored Threads My Blog My Photo Gallery My Profile My Homepage Start New Thread Enable EMail Alerts Rate this Thread
|
©2010 Advenet LLC Privacy Policy - Terms of Use
This website includes both content owned or controlled by Advenet as well as content owned or controlled by third parties.