dare_advection_diffusion.m

classdef dare_advection_diffusion < DARE.Model
% Implicit Euler discretization of a finite-difference convection-diffusion
% model.
%
% This model is constructed in such a way, that it can follow a given
% trajectory - or at least it tries to follow it. Use the SIMULATE function
% to simulate!
%
% Andreas Schmidt, 2016

    properties
        mu = [1 1 3];
        mu_names = {'mu_diffusion', 'mu_adv_y', 'lambda_tracking'};
        mu_ranges = {  [0.1, 2], [0,10], [1, 10]};
        
        n0;
        dt = 0.05;
    end
        
    methods
        function this = dare_advection_diffusion(n)
            this.n = n*n + 1;
            this.n0 = n;
            this.m = 1;
            this.p = 1;
            
            this.RB_orthonormalize_E = false;
            this.RB_error_indicator = 'normalized_residual';
            this.RB_pod_tolerance = 0.85;
            this.RB_greedy_tolerance = 1e-3;
            this.RB_M_train = ParameterSampling.Uniform(6);
        end
        
        function r = E_comp(this)
            r = {speye(this.n), ...
                 blkdiag(this.my_2d_matrix(this.n0, '0', '0', '0', '1'),0) ...
                 blkdiag(this.my_2d_matrix(this.n0, '0', '-1', '0', '0'),0)};
        end
        function r = E_coeff(this)
            r = [1, -this.dt*this.mu(1), -this.dt*this.mu(2)];
        end
        
        function r = A_comp(this)
            C = this.C();
            r = {[speye(this.n-1), zeros(this.n-1,1); -C 1]};
        end
        function r= A_coeff(this)
            r = 1;
        end
        
        function r = B_comp(this)
            v = this.my_2d_vector(this.n0, '( (x>=0.2)&(x<=0.6)&(y>=0.4)&(y<=0.6) )');
            r = {[10*1/numel(nonzeros(v))*v;0]};
        end
        function r = B_coeff(this)
            r = this.dt;
        end
        
        function r = C(this)
            v = this.my_2d_vector(this.n0, '1')';
            r = v/numel(nonzeros(v));
        end
        
        function r = C_comp(this)
            r = {[this.C(), 0], [zeros(1,this.n-1) 1]};
        end
        function r = C_coeff(this)
            r = [1 this.mu(3)];
        end
        function r = R_coeff(this)
            r = 0.01;
        end
        
        function s = get_ss(this, model_data, dsim)
            % GET_SS Get the state space model for simulation purpose
            % If dsim is provided, a closed-loop simulation is performed.
            [E,A,B,C,~,R] = this.assemble(model_data);
            outputName = {'y(t)', 'u(t)'};
            C = this.C();
            
            if exist('dsim','var') && ~isempty(dsim)
                K = inv(R+B'*dsim.Z*dsim.Z'*B)*B'*dsim.Z*dsim.Z'*A;
                A = A - B*K;
                C = [C 0; K];
            else
                C = [C 0; 0*C];
            end

            s = dss(A, [B, [zeros(this.n-1,1); 1]], C, [], full(E),this.dt, ...
                'InputName', {'u(t)', 'r(t)'}, ...
                'OutputName', outputName);
        end
        
        function [t,y,u,r,x] = simulate(this, model_data, T, r, sim, x0)
            % SIMULATE - Simulate the system for a given time
            %    T ...... final time
            %    r ...... reference value or vector
            %    sim .... simulation used for feedback
            %
            %  This function returns the following vectors:
            %    t ....... time span for simulation
            %    y ....... output value
            %    u ....... control signal
            %    r ....... the reference signal
            %    x .......
            if nargin < 5, sim = []; end
            if nargin < 6, x0 = zeros(this.n,1); end
            
            t = (0:this.dt():T)';
            if isa(r, 'function_handle'), r = r(t); end
            if length(r) == 1, r = ones(length(t),1)*r; end
            
            ss = this.get_ss(model_data, sim);
            
            [y_,t,x] = lsim(ss, [0*t, r], t, x0);
            y = y_(:,1);
            u = y_(:,2);
        end

    end
    
    methods(Access=protected)
        function [A, name] = my_2d_matrix(this,n0,fx_str,fy_str,g_str,mu_str)
        %
        %  Generates the stiffness matrix A for the finite difference
        %  discretization (equidistant grid) of the PDE
        %
        %   mu* laplace(u) - fx du/dx - fy du/dy - g u   =   r.h.s. on Omega
        %                 
        %                                        u   =   0          on dOmega
        % 
        %  Omega = (0,1)x(0,1)   (unit square). 
        %
        %  This function is just used as an easy way to generate test problems
        %  rather than to solve PDEs.
        %
        %  Calling sequence:
        %   
        %    [A, name] = fdm_2d_matrix( n0, fx_str, fy_str, g_str )
        %
        %  Input:
        %   
        %    n0        number of inner grid points in each dimension;
        %    fx_str    string describing the function fx in the space variables
        %              'x' and 'y', e.g., fx_str = 'sin(x+2*y)+3';              
        %    fy_str    string describing the function fy in the space variables
        %              'x' and 'y';              
        %    g_str     string describing the function g in the space variables
        %              'x' and 'y'. 
        %
        %  Output:
        %
        %    A         n-x-n sparse stiffness matrix, where n = n0^2;
        %    name      string describing the problem.
        %
        %
        %  LYAPACK 1.0 (Thilo Penzl, May 1999)

            % Input data not completely checked!

            na = nargin;

            if na~=5
            mu_str = '1';
            end

            name = ['FDM-2D: fx=',fx_str,'; fy=',fy_str,'; g=',g_str];

            n2 = n0*n0;

            h = 1.0/(n0+1);                            

            h2 = h*h;

            t1 = 4.0/h2;
            t2 = -1.0/h2;
            t3 = 1.0/(2.0*h);

            len = 5*n2-4*n0;
            I = zeros(len,1);
            J = zeros(len,1);
            S = zeros(len,1);
            ptr = 0;                                  % Pointer
            i = 0;                                    % Row Number

            for iy = 1:n0
            y = iy*h;
            for ix = 1:n0
            x = ix*h;

            i = i+1;
            fxv = eval(fx_str);
            fyv = eval(fy_str);
            gv = eval(g_str);
            mu = eval(mu_str);

            if iy>1
            ptr = ptr+1;                        % A(i,i-n)
            I(ptr) = i;
            J(ptr) = i-n0;
            S(ptr) = mu*t2-fyv*t3;
            end

            if ix>1
            ptr = ptr+1;                        % A(i,i-1)
            I(ptr) = i;
            J(ptr) = i-1;
            S(ptr) = mu*t2-fxv*t3;
            end 

            ptr = ptr+1;                          % A(i,i)
            I(ptr) = i;
            J(ptr) = i;
            S(ptr) = mu*t1+gv;

            if ix<n0
            ptr = ptr+1;                        % A(i,i+1)
            I(ptr) = i;
            J(ptr) = i+1;
            S(ptr) = mu*t2+fxv*t3;
            end

            if iy<n0        
            ptr = ptr+1;                        % A(i,i+n0)
            I(ptr) = i;
            J(ptr) = i+n0;
            S(ptr) = mu*t2+fyv*t3;
            end

            end  
            end

            A = -sparse(I,J,S,n2,n2);
        end

        function v = my_2d_vector(this,n0,f_str)
        %
        %  Generates a vector v which contains the values of a function f(x,y) 
        %  on an equidistant grid in the interior of the unit square. The grid
        %  points are numbered consistently with those used in the function
        %  'fdm_2d_matrix'.
        %
        %  This function is just used as an easy way to generate test problems
        %  rather than to solve PDEs.
        %
        %  Calling sequence:
        %   
        %    v = fdm_2d_vector( n0, f_str)
        %
        %  Input:
        %   
        %    n0        number of inner grid points in each dimension;
        %    f_str     string describing the function f in the space variables 'x'
        %              and 'y', e.g., f_str = 'sin(x+2*y)+3'. 
        %
        %  Output:
        %
        %    v         vector of order n = n0^2 containing the values of f(x,y).
        %
        %
        %  LYAPACK 1.0 (Thilo Penzl, May 1999)

        % Input data not completely checked!

            na = nargin;

            if na~=2
            error('Wrong number of input parameters.');
            end

            h = 1.0/(n0+1);                             

            n2 = n0*n0;

            v = zeros(n2,1);

            i = 0;

            for iy = 1:n0
            y = iy*h;
            for ix = 1:n0
            x = ix*h;
            i = i+1;
            v(i) = eval(f_str);
            end
            end  
        end
    end
end