Step function in RL environment cannot call helper functions defined below. Where is my fault?

Hey there everybody,

I'm currently working on a RL environment for my bachelor's thesis. I've started defining the environment and I'm continuously validating it. When it comes to calling helper functions ('getResistance' and 'ecms'), Matlab cannot find/call these functions that I defined in the optional methods section below.

Outputs are either "Cannot find an exact (case-sensitive) match for 'ecms'" or "Undefined function 'getResistance' for input arguments of type 'double'". Please note that I already exchanged the getResistance function call for the lines of code that are in the function's body.

The code below might help answering my question.

I am very thankful for an answer!

Max

classdef Environment < rl.env.MATLABEnvironment
    %ENVIRONMENT: Template for defining custom environment in MATLAB.
    
    %% Properties (set properties' attributes accordingly)
    properties
        dr1 = load('dr1');
        dr2 = load('dr2');

        t = 1;
        cycle = zeros(6,10);
        soc = rand(1,10);
        
        % gear ratios
        i_gear = [5.354, 3.243, 2.252, 1.636, 1.211, 1.000, 0.865, 0.717, 0.601]';
        % final drive ratio
        i_final = 3.07;
        % mass factors
        k_m = [1.32, 1.16, 1.11, 1.08, 1.07, 1.07, 1.07, 1.06, 1.06]';
        % [W] e-motor power
        P_e = 15000;
        % rolling resistance coefficient
        f_roll = 0.011;
        % [kg] vehicle mass
        m = 1625;
        % [kg] additional mass
        m_a = 0;
        % [m/s^2] gravitational acceleration
        g = 9.81;
        % [kg/m^3] air density at 20°C
        rho_air = 1.2;
        % [m^2] air drag coefficient multiplied by face area
        cda = 0.54;
        % [m] wheel radius
        r_wheel = 0.32;
        % [W] Power of attachements
        P_att = 800;
        % drivetrain efficiency
        eff_dt = 0.9;
        % [J/kg] fuel lower heat value
        H_f = 42.5e+6;
        % [kg/J] reciprocal fuel lower heat value
        e_ICE = 2.35294e-8;
        
        % battery efficiency
        eff_bat = 0.95;
        % SOC list
        soc_list = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0];
        % [As] battery capacity
        Q_b = 20*3600;
        % [Ohm] internal resistance of battery
        R_int = 0.334 / 20;
        % [V] open circuit voltage at SOC from SOC list
        V_oc  = [37.78, 45.30, 46.01, 46.56, 46.94, 47.38, 48.00, 48.78, 49.68, 50.66, 51.74];
        % reference SOC
        SOC_ref = 0.55;
        % minimal SOC
        SOC_min = 0.3;
        % maximal SOC
        SOC_max = 0.8;
        
        % [rad/s] ICE speed list (0.10472=2pi/60)
        w_ICE_list = 0.10472*[1100, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6300];
        % [Nm] ICE torque list
        t_ICE_list = [0, 30, 50, 70, 90, 110, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310, 330, 350];
        % ICE efficiency map at w_ICE and t_ICE
        eff_ICE_table = 0.01*[...
            0.01 21.18 26.06 28.71 30.80 32.33 33.35 34.16 34.57 34.16 33.22 32.70 32.21 31.96 31.96 31.96 31.96 31.96
            0.01 21.18 26.47 29.72 31.84 33.09 34.16 34.86 35.29 35.00 35.00 34.72 34.02 33.22 32.96 32.33 32.21 32.09
            0.01 21.18 26.47 29.72 31.73 32.96 34.02 34.86 35.29 35.29 35.29 35.29 35.15 34.72 34.16 33.48 33.22 33.09
            0.01 21.18 26.47 29.41 31.61 32.83 34.02 34.86 35.29 35.29 35.29 35.29 35.29 35.00 34.57 34.29 34.02 33.48
            0.01 20.17 25.67 28.42 30.91 32.33 33.48 34.29 35.00 35.29 35.29 35.29 35.29 34.86 34.29 33.88 32.58 32.21
            0.01 20.17 25.29 28.24 30.36 31.84 32.83 33.75 34.29 34.72 34.86 34.72 34.43 33.88 33.22 32.33 31.37 30.25
            0.01 20.17 24.20 27.50 29.41 30.80 31.73 32.70 33.35 33.88 33.88 33.35 32.58 32.58 31.84 30.80 30.04 29.72
            0.01 18.82 23.53 26.47 28.52 29.72 30.80 31.61 31.73 31.61 31.26 31.37 31.14 31.26 30.47 30.04 29.72 29.72
            0.01 18.82 23.21 25.51 27.32 28.52 29.72 30.25 30.36 30.25 30.04 30.04 29.62 29.51 29.41 29.41 29.41 29.41
            0.01 18.82 21.44 24.20 26.22 27.50 28.62 28.62 28.52 28.52 28.52 28.42 28.52 28.71 28.71 28.71 28.71 28.71
            0.01 18.82 20.17 23.02 24.62 25.67 26.31 26.64 26.64 25.82 25.51 24.55 24.55 24.55 24.55 24.55 24.55 24.55
            0.01 18.82 20.17 22.00 23.33 24.20 24.91 25.29 25.67 24.91 24.20 24.20 24.20 24.20 24.20 24.20 24.20 24.20];
        
        % [kg/s] consumption rate map at w_ICE and t_ICE
        consrate_ICE_table = cell2mat(struct2cell(load('consrate_ICE_table')));
        
        % [rad/s] ICE drag speed list
        w_ICE_drag_list  = 0.10472*[0, 1000, 2000, 3000, 4000, 5000, 6000];
        % [Nm] ICE drag torque list
        t_ICE_drag_list  = [6.3, 6.3, 8.4, 10.7, 13.1, 17.0, 20.7];
        
        % [rad/s] e-motor speed list
        w_em_list = 0.10472*[0, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000];
        % [Nm] e-motor torque list
        t_em_list = 15000 / 40000 *[-260, -230, -200, -170, -140, -110, -80, -50, -20, -11, -5, 0, 5, 10 ,20 ,50, 80, 110, 140, 170, 200, 230, 260];
        
        % e-motor efficiency map at w_em and t_em
        eff_em_table = 0.01*[...
            0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0 50.0
            0.01 10.0 25.0 40.0 62.0 73.6 80.0 85.0 85.0 80.0 70.0 50.0 70.0 80.0 85.0 86.0 82.5 78.7 74.6 71.6 68.8 66.4 64.8
            70.0 73.2 77.0 80.0 83.5 86.5 89.4 91.5 90.0 85.0 70.0 50.0 70.0 82.0 90.0 91.6 90.0 87.5 85.5 83.1 81.1 78.8 76.4
            80.4 82.2 84.5 86.4 88.3 90.3 91.7 93.0 90.0 83.0 70.0 50.0 70.0 82.0 90.0 92.9 92.3 91.2 89.5 87.5 86.3 85.0 84.0
            84.5 86.1 87.8 89.5 91.1 92.7 94.0 94.2 92.0 87.0 70.0 50.0 70.0 85.0 91.2 94.0 93.1 92.2 90.9 59.5 88.0 88.0 88.0
            89.5 89.5 89.5 91.0 92.4 93.8 94.5 94.7 93.0 89.0 70.0 50.0 70.0 85.0 91.0 93.4 93.1 92.2 91.0 90.0 90.0 90.0 90.0
            91.5 91.5 91.5 91.5 92.8 94.1 94.8 95.3 93.0 88.0 70.0 50.0 70.0 84.0 90.5 93.1 92.8 92.1 91.0 91.0 91.0 91.0 91.0
            92.8 92.8 92.8 92.8 92.8 94.1 95.0 95.7 92.7 87.0 70.0 50.0 70.0 82.0 90.0 92.7 92.4 91.8 91.8 91.8 91.8 91.8 91.8
            92.6 92.6 92.6 92.6 92.6 93.8 94.7 95.3 92.0 85.0 70.0 50.0 70.0 81.0 89.0 92.2 92.0 91.0 91.0 91.0 91.0 91.0 91.0
            93.2 93.2 93.2 93.2 93.2 93.2 94.3 94.7 92.0 85.0 70.0 50.0 70.0 80.0 87.5 91.6 91.0 91.0 91.0 91.0 91.0 91.0 91.0
            93.0 93.0 93.0 93.0 93.0 93.0 94.0 94.4 90.4 82.0 70.0 50.0 70.0 80.0 86.5 91.0 90.5 90.5 90.5 90.5 90.5 90.5 90.5
            93.5 93.5 93.5 93.5 93.5 93.5 93.5 93.8 90.0 82.0 70.0 50.0 70.0 78.0 85.6 90.3 90.0 90.0 90.0 90.0 90.0 90.0 90.0
            92.8 92.8 92.8 92.8 92.8 92.8 92.8 93.0 90.0 80.0 70.0 50.0 70.0 75.0 85.0 88.0 88.0 88.0 88.0 88.0 88.0 88.0 88.0
            ];
    end
    
    properties
        % Initialize system state [v, a, alpha, gear, soc]'
        State = zeros(5,1)
    end
    
    properties(Access = protected)
        % Initialize internal flag to indicate episode termination
        IsDone = false
    end
    
    %% Necessary Methods
    methods
        % Contructor method creates an instance of the environment
        % Change class name and constructor name accordingly
        function this = Environment()
            % Initialize Observation settings
            ObservationInfo = rlNumericSpec([5 1]);
            ObservationInfo.Name = 'Vehicle State';
            ObservationInfo.Description = 'v, a, alpha, gear, soc';
            
            % Initialize Action settings
            ActionInfo = rlNumericSpec([2 1]);
            ActionInfo.Name = 'Vehicle Action';
            ActionInfo.Description = 'lambda, sigma';
            
            % The following line implements built-in functions of RL env
            this = [e-mail](ObservationInfo,ActionInfo);
            
            % Initialize property values and pre-compute necessary values
            % updateActionInfo(this);
        end
        
        % Apply system dynamics and simulates the environment with the
        % given action for one step.
        function [Observation,Reward,IsDone,LoggedSignals] = step(this,Action)
            LoggedSignals = [];
            
            this.soc(1) = this.SOC_ref;
            
            v = this.State(1);
            a = this.State(2);
            alpha = this.State(3);
            gear = this.State(4);
            soc = this.State(5);
            
            % calculate driving resistance
            F_air   = 0.5 * this.rho_air * this.cda * v^2; % [N] air resistance
                F_roll  = this.f_roll * (this.m + this.m_a) * this.g * cos(alpha); % [N] rolling resistance
                F_incl  = (this.m + this.m_a) * this.g * sin(alpha); % [N] incline resistance
                F_acc   = (this.k_m(gear) * this.m + this.m_a) * a; % [N] required force for acceleration
            
            F_w     = F_air + F_roll + F_incl + F_acc; % [N] total driving resistance
            
            w_wheel = v / this.r_wheel; % [rad/s] wheel speed
            
            w_crank = w_wheel * this.i_gear(gear) * this.i_final; % [rad/s] crank speed
            
                t_wheel = F_w * this.r_wheel; % [Nm] wheel torque
                P_wheel = w_wheel * t_wheel;
                P_crank = (P_wheel > 0) * P_wheel / this.eff_dt...
                        + (P_wheel == 0) * 0 + ...
                        + (P_wheel < 0) * P_wheel * this.eff_dt;
            if w_crank > 0
            t_crank = P_crank / w_crank; % [Nm] crank torque
            else
            t_crank = 0; % [Nm] crank torque
            end
            
            [t_em_opt, fuel_opt, w_ICE, t_ICE] = ecms(t_crank, w_crank, Action(1), Action(2), soc);
            
            if this.t < length(this.cycle.v)
                new_t = this.t + 1;
                this.t = new_t;
                IsDone = false;
            elseif this.t == length(this.cycle.v)
                IsDone = true;
            end
            
            v = this.cycle.v(this.t);
            a = this.cycle.a(this.t);
            alpha = this.cycle.alpha(this.t);
            gear = this.cycle.gear(this.t);
            this.soc(this.t) = soc ;
            
            Observation = [v; a; alpha; gear; soc];
            this.State = Observation;
        end
        
        % Reset environment to initial state and output initial observation
        function [InitialObservation, t, cycle] = reset(this)
            this.t = 1;
            this.cycle = this.dr1.dr1.dr1_6;
            
            v = this.cycle.v(1);
            a = this.cycle.a(1);
            alpha = this.cycle.alpha(1);
            gear = this.cycle.gear(1);
            soc = this.SOC_ref;
            
            InitialObservation = [v; a; alpha; gear; soc];
            this.State = InitialObservation;
            this.t = 1;
            
            % (optional) use notifyEnvUpdated to signal that the
            % environment has been updated (e.g. to update visualization)
            % notifyEnvUpdated(this);
        end
    end
    %% Optional Methods (set methods' attributes accordingly)
    methods (Access = public)
        function [F_w, w_wheel, w_crank, t_crank] = getResistance(v, a, alpha, gear)
                F_air   = 0.5 * this.rho_air * this.cda * v^2; % [N] air resistance
                F_roll  = this.f_roll * (this.m + this.m_a) * this.g * cos(alpha); % [N] rolling resistance
                F_incl  = (this.m + this.m_a) * this.g * sin(alpha); % [N] incline resistance
                F_acc   = (this.k_m(gear) * this.m + this.m_a) * a; % [N] required force for acceleration
            
            F_w     = F_air + F_roll + F_incl + F_acc; % [N] total driving resistance
            
            w_wheel = v / this.r_wheel; % [rad/s] wheel speed
            
            w_crank = w_wheel * this.i_gear(gear) * this.i_final; % [rad/s] crank speed
            
                t_wheel = F_w * this.r_wheel; % [Nm] wheel torque
                P_wheel = w_wheel * t_wheel;
                P_crank = (P_wheel > 0) * P_wheel / this.eff_dt...
                        + (P_wheel == 0) * 0 + ...
                        + (P_wheel < 0) * P_wheel * this.eff_dt;
            if w_crank > 0
            t_crank = P_crank / w_crank; % [Nm] crank torque
            else
            t_crank = 0; % [Nm] crank torque
            end
        end
        
        function [t_em_opt, fuel_opt, w_ICE, t_ICE] = ecms(t_crank, w_crank, lambda, sigma, soc)
            dt = 1; % [Nm] torque step
            lambda = lambda; % equivalence factor
            
            % [Nm] max e-motor torque (powering mode)
            T_e_maxp = min(260, this.P_e/w_crank); 
            % [Nm] max e-motor torque (generating mode)
            T_e_maxg = min(260, 1.25*this.P_e/w_crank); 
            % [Nm] ICE drag torque
            t_ICE_drag = interp1(this.w_ICE_drag_list,this.t_ICE_drag_list,w_crank,'linear','extrap');
            
            % [Nm] max ICE torque
            if w_crank < 0.10472 * 1250
                t_ICE_max = min(350-5.41*(0.10472*1250-w_crank), t_crank + T_e_maxg);
            elseif w_crank >= 0.10472 * 1250 && w_crank < 0.10472 * 4000
                t_ICE_max = min(350, t_crank + T_e_maxg);
            else
                t_ICE_max = min(350-0.54*(w_crank-0.10472*4000), t_crank + T_e_maxg);
            end
            
            % [Nm] min ICE torque
            if t_crank < T_e_maxp - t_ICE_drag
                t_ICE_min = 0;
            elseif t_crank < T_e_maxp
                t_ICE_min = 1;
            else
                t_ICE_min = t_crank - T_e_maxp;
            end
            
            % initialization: set first falue of fuel_opt to large number.
            fuel_opt = 1; % [kg/s]
            fuel_min = 1; % [kg/s]
            t_em_opt = 0; % [Nm]
            
            % modify engine speed to delivers power at min. 1100 rpm
            % in case that w_crank < 1100 rpm, the clutch slips
            w_ICE = max(0.10472*1100, w_crank);
            
            %%% ECMS algorithm %%%
            if w_crank == 0
                fuel_opt = 0;
                t_em_opt = 0;
                
            elseif w_crank > 0 && w_crank < 115.192 && t_crank > 0 % vehicle launches
                % slipping clutch launch (hybrid launch possible)
                for t_ICE = t_ICE_min : dt : t_ICE_max % search for optimal ICE/EM ratio
                    
                    if t_ICE == 0
                        t_em = t_crank + t_ICE_drag; % [Nm]
                    else
                        t_em = t_crank - t_ICE; % [Nm]
                    end
                    
                    % [kg/s] current consumption rate
                    fuel_h = this.e_ICE * 115.192 * t_ICE / interp2(this.t_ICE_list, this.w_ICE_list, this.eff_ICE_table, t_ICE, 115.192) + ... 
                        lambda * (t_em>=0) * 115.192 * t_em / interp2(this.t_em_list, this.w_em_list, this.eff_em_table, t_em, 115.192) + ...
                        lambda * (t_em <0) * 115.192 * t_em * interp2(this.t_em_list, this.w_em_list, this.eff_em_table, t_em, 115.192);
                    if fuel_h < fuel_min
                        fuel_min = fuel_h;
                        % [kg/s] optimal current consumption rate
                        fuel_opt = this.e_ICE * 115.192 * t_ICE / interp2(this.t_ICE_list, this.w_ICE_list, this.eff_ICE_table, t_ICE, 115.192);
                        % [Nm] optimal e-motor torque
                        t_em_opt = t_em;
                    end
                end
                
                % check feasibility of closed clutch launch (electric only launch)
                t_em_l = t_crank + t_ICE_drag;
                if t_ICE_min == 0 && t_em_l <= T_e_maxp
                    fuel_e = lambda * w_crank * t_em_l / interp2(this.t_em_list, this.w_em_list, this.eff_em_table, t_em_l, w_crank);
                    if fuel_e < fuel_min
                        fuel_opt = 0;
                        t_em_opt = t_em_l;
                    end
                end
                
            elseif w_crank > 0 && w_crank < 115.192 && t_crank <= 0 && t_crank > - T_e_maxg - t_ICE_drag %%% vehicle stops
                % electric charge only
                fuel_opt = 0;
                t_em_opt = t_crank + t_ICE_drag;
                
            elseif w_crank > 0 && w_crank < 115.192 && t_crank <= - T_e_maxg - t_ICE_drag
                % electric charge only
                fuel_opt = 0;
                t_em_opt = - T_e_maxg;
                
            elseif w_crank >= 115.192 && t_crank <= - T_e_maxg
                % electric charge only
                fuel_opt = 0;
                t_em_opt = max((t_crank + t_ICE_drag),- T_e_maxg);
                
            elseif w_crank >= 115.192 && t_crank > - T_e_maxg
                % hybrid mode possible
                for t_ICE = t_ICE_min : dt : t_ICE_max
                    
                    if t_ICE == 0
                        t_em = t_crank + t_ICE_drag;
                    else
                        t_em = t_crank - t_ICE;
                    end
                    
                    fuel_h = this.e_ICE * w_crank * t_ICE / interp2(this.t_ICE_list, this.w_ICE_list, this.eff_ICE_table, t_ICE, w_crank) + ...
                        ef_RL * (t_em>=0) * w_crank * t_em / interp2(this.t_em_list, this.w_em_list, this.eff_em_table, t_em, w_crank) + ...
                        ef_RL * (t_em <0) * w_crank * t_em * interp2(this.t_em_list, this.w_em_list, this.eff_em_table, t_em, w_crank);
                    if fuel_h < fuel_min
                        fuel_min = fuel_h;
                        fuel_opt = this.e_ICE * w_ICE * t_ICE / interp2(this.t_ICE_list, this.w_ICE_list, this.eff_ICE_table, t_ICE, w_ICE);
                        t_em_opt = t_em;
                    end
                end
            end
            
            t_ICE = t_crank - t_em_opt;
        end

        function reward = getreward(sigma, w_ICE, t_ICE, delta_ess, last_ess, deltatorque, DeltaTorque, soc, P_batt_10)
                % [kg] last second's fuel consumption of running engine
                cons_run        = sigma * w_ICE * t_ICE / (this.H_f * interp2(this.t_ICE_list, this.w_ICE_list, this.eff_ICE_table, t_ICE, w_ICE)); 
                % [kg] fuel consumption for engine start (0.3g per start)
                cons_start      = (delta_ess==1) * 0.0003 * delta_ess; 
                % penalty for frequent engine starts/stops
                startstop       = (last_ess<=15) * (0.5 * abs(delta_ess) * cos(pi/15 * last_ess) + 0.5); 
                % penalty for large engine torque jumps (turbocharged engines)
                icetorque       = max(0, deltatorque/DeltaTorque - 0.5); 
                % penalty for difference from reference SOC (0.55)
                battery_soc     = (soc-0.55)^2/((0.8-soc)*(soc-0.3)); 
                % penalty for large battery power
                battery_power   = max(0, rms(P_batt_10) - P_batt_max/2);
            % total reward
                C1 = 1;
                C2 = 2;
                C3 = 4;
                C4 = 2/P_batt_max;
            reward = cons_run + cons_start + C1*startstop + C2*icetorque + C3*battery_soc + C4*battery_power;
        end
    end
    
    methods (Access = protected)
        % (optional) update visualization everytime the environment is updated
        % (notifyEnvUpdated is called)
        function envUpdatedCallback(this)
        end
    end
end