ODE-Test-Problems/src/+otp/+pendulumdae/PendulumDAEProblem.m at 2d282cceb1a525e0919c0f67a0ca420e1a9efd52 · ComputationalScienceLaboratory/ODE-Test-Problems · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
classdef PendulumDAEProblem < otp.Problem
    %CONSTRAINEDPENDULUM PROBLEM This is a Hessenberg Index-2 DAE problem posed in terms of three constraints
    %

    properties (SetAccess=protected)
        RHSDifferential
        RHSAlgebraic
        RHSHalfExplicit
    end

    properties (Dependent)
      Y0Differential
      Y0Algebraic
    end

    methods
        function obj = PendulumDAEProblem(timeSpan, y0, parameters)
            obj@otp.Problem('Constrained Pendulum', 7, timeSpan, y0, parameters);
        end
    end

    methods (Access = protected)
        function onSettingsChanged(obj)
            m = obj.Parameters.Mass;
            l = obj.Parameters.Length;
            g = obj.Parameters.Gravity;

            % get initial energy
            initialconstraints = otp.pendulumdae.fAlgebraic([], obj.Y0Differential, g, m, l, 0);
            E0 = initialconstraints(3);

            % 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', otp.pendulumdae.mass([], [], g, m, l, E0), ...
                'Jacobian', @(t, y) otp.pendulumdae.jacobian(t, y, g, m, l, E0), ...
                'JacobianVectorProduct', ...
                @(t, y, v) otp.pendulumdae.jacobianVectorProduct(t, y, v, g, m, l, E0), ...
                'Vectorized', 'on');

            % Generate the constituent RHS for the differential part
            obj.RHSDifferential = otp.RHS(@(t, y) otp.pendulumdae.fDifferential(t, y, g, m, l, E0), ...
                'Jacobian', @(t, y) otp.pendulumdae.jacobianDifferential(t, y, g, m, l, E0), ...
                'JacobianVectorProduct', ...
                @(t, y, v) otp.pendulumdae.jacobianVectorProductDifferential(t, y, v, g, m, l, E0), ...
                'Vectorized', 'on');

            % Generate the constituent RHS for the algebraic part
            obj.RHSAlgebraic = otp.RHS(@(t, y) otp.pendulumdae.fAlgebraic(t, y, g, m, l, E0), ...
                'Jacobian', @(t, y) otp.pendulumdae.jacobianAlgebraic(t, y, g, m, l, E0), ...
                'JacobianVectorProduct', ...
                @(t, y, v) otp.pendulumdae.jacobianVectorProductAlgebraic(t, y, v, g, m, l, E0), ...
                'JacobianAdjointVectorProduct', ...
                @(t, y, v) otp.pendulumdae.jacobianAdjointVectorProductAlgebraic(t, y, v, g, m, l, E0, v), ...
                'HessianAdjointVectorProduct', ...
                @(t, y, v, u) otp.pendulumdae.hessianAdjointVectorProductAlgebraic(t, y, v, u, g, m, l, E0), ...
                'Vectorized', 'on');

            % Generate an Auxillary RHS for the half-explicit Jacobian
            obj.RHSHalfExplicit = otp.RHS([], ...
                'Jacobian', @(t, y) otp.pendulumdae.jacobianHalfExplicit(t, y, g, m, l, E0), ...
                'JacobianVectorProduct', ...
                @(t, y, v) otp.pendulumdae.jacobianVectorProductHalfExplicit(t, y, v, g, m, l, E0), ...
                'Vectorized', 'on');

        end

        function label = internalIndex2label(~, index)
            switch index
                case 1
                    label = 'x position';
                case 2
                    label = 'y position';
                case 3
                    label = 'x velocity';
                case 4
                    label = 'y velocity';
                case 5
                    label = 'Control 1';
                case 6
                    label = 'Control 2';
                case 7
                    label = 'Control 3';
            end
        end

        function sol = internalSolve(obj, varargin)
            sol = internalSolve@otp.Problem(obj, 'Solver', @otp.utils.solvers.dae34, varargin{:});
        end
    end

    methods
        function y0differential = get.Y0Differential(obj)
            y0differential = obj.Y0(1:4);
        end

        function y0algebraic = get.Y0Algebraic(obj)
            y0algebraic = obj.Y0(5:end);
        end
    end

end