Search code, repositories, users, issues, pull requests...

Copy link

Member

@AndreyAPopov AndreyAPopov commented May 10, 2022

For use in half-explicit methods and for data assimilation.


 Initial Commit of the constrained pendulum problem

9c47f78

For use in half-explicit methods and for data assimilation.

@AndreyAPopov AndreyAPopov requested review from Steven-Roberts, amitns122 and reid-g

May 10, 2022 15:56

AndreyAPopov added 2 commits

May 10, 2022 14:44


 very useful fix for lambda

6d54390


 added Z0 and fixed jacobia

0caf7cc

Steven-Roberts reviewed

src/+otp/+constrainedpendulum/+presets/Canonical.m Outdated

Comment on lines 5 to 6

% Ascher, Uri. "On symmetric schemes and differential-algebraic equations."

% SIAM journal on scientific and statistical computing 10.5 (1989): 937-949.

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this the right source? I don't see anything on a pendulum.

Copy link

Member

@amitns122 amitns122 May 11, 2022 •

edited

Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hairer, E., Roche, M., Lubich, C. (1989). Description of differential-algebraic problems. In: The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods. Lecture Notes in Mathematics, vol 1409. Springer, Berlin, Heidelberg. https://doi.org/10.1007/BFb0093948

That is the pendulum problem source.

The Ascher book has a nice formulation as a Hessenberg indeex 2 DAE( linear algebraic variables) we use.

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah the citation is wrong.

Steven-Roberts reviewed

src/+otp/+constrainedpendulum/ConstrainedPendulumParameters.m Outdated

%Gravity defines the magnitude of acceleration towards the ground

Gravity %MATLAB ONLY: (1,1) {mustBeNonnegative}

%Mass defines the mass of the pendulum

Mass %MATLAB ONLY: (1,1) {mustBeNonnegative}

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sharper validation: {mustBeReal, mustBeFinite, mustBePositive}

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Better yet, it should match https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/3d4cb2262931eaf261ca3061479b81928b1ef419/src/+otp/+pendulum/PendulumParameters.m

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea

Steven-Roberts reviewed

src/+otp/+constrainedpendulum/ConstrainedPendulumProblem.m Outdated

@@ -0,0 +1,35 @@

classdef ConstrainedPendulumProblem < otp.Problem

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps just PendulumDAEProblem.

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This class ought to override internalIndex2label

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree with both

Steven-Roberts reviewed

src/+otp/+constrainedpendulum/ConstrainedPendulumProblem.m Outdated

%CONSTRAINEDPENDULUM PROBLEM This is a Hessenberg Index-2 DAE problem posed in terms of three constraints

properties

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These should have SetAccess = private or protected

Steven-Roberts reviewed

src/+otp/+constrainedpendulum/ConstrainedPendulumProblem.m Outdated

initialconstraints = otp.constrainedpendulum.constraints([], obj.Y0, g, m, l, 0);

E0 = initialconstraints(3);

obj.RHS = otp.RHS(@(t, y) otp.constrainedpendulum.f(t, y, g, m, l, E0));

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This RHS will lead to confusion as a call to solve will not enforce constraints, right? Ideally, RHS should be in a form that a built-in solver can integrate, but we can let that slide for an index-2 DAE. I suggest the default RHS should use a mass matrix and an f that operates on the full state (diff and alg). Then there can be a RHSDifferential and RHSAlgebraic to get the individual pieces.

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also it looks like 'Vectorized', 'on' can be added to RHS

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with this

Steven-Roberts reviewed

src/+otp/+constrainedpendulum/constraintsjacobian.m Outdated

@@ -0,0 +1,27 @@

function dc = constraintsjacobian(~, state, g, m, ~, ~)

Copy link

Member

@Steven-Roberts Steven-Roberts May 11, 2022 •

edited

Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What Jacobian is this exactly? For the DAE
y' = f(y, z)
0 = g(y)
is it g_y or f_z*g_y

Copy link

Member

@amitns122 amitns122 May 11, 2022 •

edited

Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

in this formulation. y' = fhat(y) - ((g_y)^T)*z. T is transpose.

so g_y * f_z = g_y*((g_y)^T)

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this jacobian is actually just g_y

@Steven-Roberts

Copy link

Member

Steven-Roberts commented May 11, 2022

Long term, I would like to have a problem class for constrained mechanical systems (Solving ODEs II pg. 464), and the pendulum would simply be a preset. This would also unify the handling of reduced index forms.

Copy link

Member Author

AndreyAPopov commented May 11, 2022

Long term, I would like to have a problem class for constrained mechanical systems (Solving ODEs II pg. 464), and the pendulum would simply be a preset. This would also unify the handling of reduced index forms.

I agree longterm. Right now we need this problem for something. We might also add a constrained shallow water (2D)

AndreyAPopov added 3 commits

May 11, 2022 09:45


 file renaming

afdfa7e


 Fixes to some of steven's complains

cd25494


 access fix

f072d66

reid-g

reid-g reviewed

src/+otp/+constrainedpendulum/+presets/Canonical.m Outdated

y0 = [sqrt(2)/2; sqrt(2)/2; 0; 0];

obj = obj@otp.constrainedpendulum.ConstrainedPendulumProblem(tspan, y0, params);

Copy link

Member

@reid-g reid-g May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will the solve functionality be disabled for this problem?

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is an interesting question. In theory we should be able to add our own small method to OTP to handle these cases

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point. Here's one idea: Add a HighIndexDAE property to https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/master/src/%2Botp/%2Butils/Solver.m which simply throws an error message that high-index DAEs are not supported by matlab solvers. Eventually when we have a constrained mechanical system problem, we will support RHS reduced to index-1 which can rely on built-in solvers.

src/+otp/+constrainedpendulum/constraintsjacobian.m Outdated

@@ -0,0 +1,27 @@

function dc = constraintsjacobian(~, state, g, m, ~, ~)

Copy link

Member

@reid-g reid-g May 11, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this to be the standard for otp when getting the Jacobian for Half-Explicit methods? Or is this an internal functionality of the problem?

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if the problem is posed in terms of invariants, this should be standard. I might rename it to invariantJacobian

Steven-Roberts reviewed

src/+otp/+pendulumdae/fdifferential.m Outdated

@@ -0,0 +1,19 @@

function dstate = fdifferential(~, state, g, m, ~, ~)

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The convention for functions is camel case now: https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/master/docs/CONTRIBUTING.md#functions

Steven-Roberts reviewed

ODE-Test-Problems/src/+otp/+zlakinetics/ZlaKineticsProblem.m

src/+otp/+pendulumdae/PendulumDAEProblem.m Outdated

% The right hand size in terms of x, y, x', y', and three control parameters z

obj.RHS = otp.RHS(@(t, y) otp.pendulumdae.f(t, y, g, m, l, E0), ...

'Mass', @(t, y) otp.pendulumdae.mass(t, y, g, m, l, E0));

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the mass matrix is constant, it can be a matrix instead of a function handle. See

Line 19 in 3d4cb22

'Mass', otp.zlakinetics.mass([], [], k, K, klA, Ks, pCO2, H), ...

for example

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was just copying the pendulum problem, which had @(t,y). I will change

AndreyAPopov added 3 commits

May 12, 2022 08:05


 final fixes

dfee988


 Added a bunch of derivatives

13cf20d


 added vectorization

0bc36d2

Copy link

Member Author

AndreyAPopov commented May 12, 2022

Unless I am missing something, I think I have addressed all the issues.

Steven-Roberts reviewed

src/+otp/+pendulumdae/PendulumDAEProblem.m Outdated

Comment on lines 11 to 12

Y0Differential

Z0Algebraic

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not opposed to having these, but if we do include this, I think all DAE problems should have these for consistency. Also, I think Y0Algebraic would be better than Z0Algebraic.

Copy link

Member

@reid-g reid-g May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't really differentiate the two. Algebraic implies there is no direct differential equation thus it is a variable obtained from differentiation of the constraint. I think it should stay.

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022 •

edited

Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this naming is pretty clear. It's mostly motivated by established naming conventions. For partitioned problems, we use <property><partition name>. E.g.

RHSStiff and RHSNonstiff
jacobianDifferential and jacobianAlgebraic
Y0Position and Y0Velocity

This is nice for autocompletion as well.

Copy link

Member

@reid-g reid-g May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My bad I misinterpreted what was being proposed. I thought you wanted Y0Differential -> Y0Algebraic not Z0Algebraic -> Y0Algebraic. I agree on the convention you proposed.

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was thinking of changing it to that already, glad to see we all think alike

Steven-Roberts reviewed

src/+otp/+pendulumdae/PendulumDAEProblem.m Outdated

% Generate the constituent RHS for the differential part

obj.RHSDifferential = otp.RHS(@(t, y) otp.pendulumdae.fDifferential(t, y, g, m, l, E0), ...

'JacobianVectorProduct', ...

@(t, y, v) otp.pendulumdae.jacobianDifferentialVectorProduct(t, y, g, m, l, E0, v), ...

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

jacobianVectorProductDifferential is the naming convention https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems/blob/master/docs/CONTRIBUTING.md#functions-for-rhs

Steven-Roberts reviewed

src/+otp/+pendulumdae/PendulumDAEProblem.m Outdated

'Vectorized', 'on');

% Generate the constituent RHS for the algebraic part

obj.RHSAlgebraic = otp.RHS(@(t, y) otp.pendulumdae.invariants(t, y, g, m, l, E0), ...

Copy link

Member

@Steven-Roberts Steven-Roberts May 12, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The terms "algebraic" and "invariant" are used inconsistently. I think "algebraic" is the more standard DAE term, so I say we should name the functions
fAlgebraic
jacobianAlgebraic
...

@Steven-Roberts

Copy link

Member

Steven-Roberts commented May 12, 2022

Unless I am missing something, I think I have addressed all the issues.

Added some final pedantic comments

AndreyAPopov added 2 commits

May 12, 2022 17:42


 name changes to make everything consistent

42d822f


 one more name fix

9b2ee98

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m Outdated

reltol = odeget(options, 'RelTol', 1e-3);

abstol = odeget(options, 'AbsTol', 1e-6);

J = odeget(options, 'Jacobian', @(t, y) jacapprox(f, t, y));

M = odeget(options, 'Mass', speye(numel(y0)));

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should determine if the mass matrix needs to be sparse from the Jacobian. Ideally, we could call eye(n, 'like', J) but this is not supported in Octave.

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m

h = tend - tc;

end

stages = zeros(numel(yc), stagenum);

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to reallocate every step

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m Outdated

np = inf;

ntol = 1e-6;

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be based on AbsTol and RelTol

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m Outdated

nmaxits = 5;

its = 0;

while norm(np) > ntol || its < nmaxits

Mc = M(staget, yc + stagedy + gh*newtonk0);

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the mass matrix is state dependent, a Newton iteration would need dM/dy. We should parse the MStateDependence property to handle this better

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m Outdated

Mc = M(staget, yc + stagedy + gh*newtonk0);

g = Mc*newtonk0 - f(staget, yc + stagedy + gh*newtonk0);

H = Mc - gh*J(staget, yc + stagedy + gh*newtonk0);

[np, ~] = gmres(H, g);

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A direct solver would be more standard. Unfortunately, Octave does not support decomposition, so lu is probably best.

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 15, 2022 •

edited

Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I found that the matrix is often very poorly scaled, so a direct solve was throwing tons of warnings.

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m Outdated

while norm(np) > ntol || its < nmaxits

Mc = M(staget, yc + stagedy + gh*newtonk0);

g = Mc*newtonk0 - f(staget, yc + stagedy + gh*newtonk0);

H = Mc - gh*J(staget, yc + stagedy + gh*newtonk0);

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022 •

edited

Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can upgrade Newton iterations to algorithm on pg 119 of Solving ODEs II. This will evaluate the Jacobian only as needed.

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m Outdated

Mc = M(tc + h, yc);

err = rms((Mc*(yc-yhat))./sc);

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rms not supported in Octave

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure why error is multiplied by Mc

Copy link

Member Author

@AndreyAPopov AndreyAPopov May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can do manual rms. Error is multiplied by Mc, as one way to do error is to simply look at the differential variables and ignore the algebraic.

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae34.m

end

function J = jacapprox(f, t, y)

Copy link

Member

@Steven-Roberts Steven-Roberts May 15, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We'll need something a bit more robust or not support finite differences at all.

AndreyAPopov added 14 commits

May 15, 2022 20:03


 Some tests.

c51c478


 moved to symmlq

46a4436


 Moved to LSQR

d23d13d


 added a bunch of stuff to the SDIRK method


 Fixed most issues with dae34, and renamed it

7ec54c6


 First commit of ESDIRK method.

c5cce23


 replace every solver with bicg

ff96669


 optimize embedded coefficients

fc51b18


 Working codes

3c01164


 Newton's method control for Lobatto

e04aa6c


 Increase the number of iterations

5a2cef0


 slight changes to dae43

73a155a


 Added Jacobian storage and recomputation criteria. Added a starting a...

026c815

...lgorithm.
The starting algorithm for FIRK can actually just use the previous stages. This speeds up computation and just works. For SDIRK this does not work, just for some Lobatto methods


 Added sparse support and support for Jacobian decompositions

0e45849

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae43.m Outdated

reltol = odeget(options, 'RelTol', 1e-3);

abstol = odeget(options, 'AbsTol', 1e-6);

J = odeget(options, 'Jacobian', []);

M = odeget(options, 'Mass', speye(numel(y0)));

Copy link

Member

@Steven-Roberts Steven-Roberts May 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Default value should be []

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae43.m Outdated

% Compute initial step size

if isempty(h)

sc = abstol + reltol*abs(y0);

Copy link

Member

@Steven-Roberts Steven-Roberts May 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

RelTol could be a vector. Should use .*

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae43.m Outdated

step = 1;

Mc = M(tspan(1), y0);

Mfull = zeros(n*stagenum, n*stagenum, 'like', Mc);

Copy link

Member

@Steven-Roberts Steven-Roberts May 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unfortunately, 'like' option is not supported in Octave

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae43.m Outdated

end

function J = jacapprox(f, t, y)

Copy link

Member

@Steven-Roberts Steven-Roberts May 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems like a useful enough to warrant it's own publicly-accessible function. We can also look into more sophisticated and robust approaches.

Steven-Roberts reviewed

src/+otp/+utils/+solvers/dae43.m

M = @(t, y) eye(numel(y0));

end

elseif ~isa(M, 'function_handle')

M = @(t, y) M;

Copy link

Member

@Steven-Roberts Steven-Roberts May 22, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We also need to support mass matrices that are only a function of time, e.g., M(t)


 Made all changes that steven requested

2d282cc

Copy link

Member Author

AndreyAPopov commented May 23, 2022

All resolved.