Initial conditions solve failed to converge. When trying to advance one time step, nonlinear solver failed to converge due to Linear Algebra error.

Hi, I'm working on a solar coupled PEM electrolyser, I imported MPPT solar charge controller model and PEM electrolysis system from Matlab add-on library. I've applied gas properties and all other variables required by MEA and all initial conditions have been set to low priority. The below is the exact error state:

Error:

An error occurred during simulation and the simulation was terminated
    Caused by:
        ['charmaine_solar_electrolyser/Solver Configuration']: Initial conditions solve failed to converge.
            When trying to advance one time step, nonlinear solver failed to converge due to Linear Algebra error.

Workspace script:

% ===================================================================
% MASTER SETUP SCRIPT (SELF-CONTAINED WITH CUSTOM UNITS)
% ===================================================================
% --- Step 1: Start with a clean slate ---
clear all;
clc;
disp('Workspace cleared. Starting master setup script...');
% --- Step 2: Define Gas Properties Directly ---
disp('Defining gas properties...');
% Data for Nitrogen (N2)
N2_T  = [-73.15; -23.15; 26.85; 76.85; 126.85; 176.85; 226.85]; % Temperature (degC)
N2_h  = [200.9; 251.0; 301.2; 351.5; 402.0; 452.7; 503.6]; % Enthalpy (kJ/kg)
N2_mu = [13.3; 16.0; 18.5; 20.8; 22.9; 24.9; 26.9]; % Viscosity (uPa*s)
N2_k  = [18.1; 22.3; 26.2; 30.0; 33.7; 37.1; 40.5]; % Conductivity (mW/(m*K))
N2_R  = 296.8; % Specific gas constant (J/kg*K)
% Data for Oxygen (O2)
O2_T  = [-73.15; -23.15; 26.85; 76.85; 126.85; 176.85; 226.85]; % Temperature (degC)
O2_h  = [184.5; 230.0; 276.0; 322.5; 369.5; 417.0; 465.0]; % Enthalpy (kJ/kg)
O2_mu = [14.8; 17.8; 20.6; 23.2; 25.6; 27.9; 30.1]; % Viscosity (uPa*s)
O2_k  = [18.3; 22.3; 26.6; 30.8; 34.9; 38.9; 42.8]; % Conductivity (mW/(m*K))
O2_R  = 259.8; % Specific gas constant (J/kg*K)
O2_D  = 20600; % Mass diffusivity (mm^2/s)
% Data for Hydrogen (H2)
H2_T  = [-73.15; -23.15; 26.85; 76.85; 126.85; 176.85; 226.85]; % Temperature (degC)
H2_h  = [2770; 3440; 4110; 4790; 5470; 6150; 6830]; % Enthalpy (kJ/kg)
H2_mu = [6.13; 7.42; 8.65; 9.81; 10.92; 11.98; 13.01]; % Viscosity (uPa*s)
H2_k  = [129; 155; 181; 205; 228; 250; 272]; % Conductivity (mW/(m*K))
H2_R  = 4124.2; % Specific gas constant (J/kg*K)
H2_D  = 70000; % Mass diffusivity (mm^2/s)
% Data for Water Vapor (H2O)
H2O_h     = [2360; 2440; 2550; 2650; 2745; 2840; 2940];   % Specific Enthalpy (kJ/kg)
H2O_mu    = [7.0; 9.0; 11.2; 13.4; 15.5; 17.6; 19.6];     % Dynamic Viscosity (uPa*s)
H2O_h_vap = [2838; 2650; 2438; 2316; 2182; 2031; 1850];   % Enthalpy of Vaporization (kJ/kg)
H2O_k     = [12.2; 15.8; 19.5; 23.5; 27.8; 32.5; 37.5];   % Thermal Conductivity (mW/(m*K))
H2O_p_sat = [1.0e-7; 7.75e-5; 0.0035; 0.0416; 0.246; 0.933; 2.64]; % Saturation Pressure (MPa)
H2O_D     = 25; % Mass diffusivity of water vapor in air (mm^2/s)
% --- Step 3: Create Solar and Load Profiles ---
solar_profile_time = (0:60:86400)';
drive_cycle_time = solar_profile_time; 
sunrise = 6 * 3600;
sunset  = 18 * 3600;
daylight_duration = sunset - sunrise;
solar_irradiance = zeros(size(solar_profile_time));
daylight_indices = solar_profile_time >= sunrise & solar_profile_time <= sunset;
daylight_time_elapsed = solar_profile_time(daylight_indices) - sunrise;
solar_irradiance(daylight_indices) = 1000 * sin(pi * daylight_time_elapsed / daylight_duration);
solar_profile = [solar_profile_time, solar_irradiance];
% Create a time vector from 0 to 10 seconds
time = (0:0.1:10)'; 
% Create an irradiance data vector (e.g., a ramp up and down)
irradiance_vals = [linspace(500,1000,51), linspace(1000,500,50)]';
% Prepend a ramp from 0 to the starting value over the first 0.01s
time = [0; 0.01; time+0.01]; % Add starting points to time
irradiance_vals = [0; 500; irradiance_vals]; % Add starting points to data
% Combine them into a two-column matrix for Simulink
irradiance_profile = [time, irradiance_vals];
% --- Step 4: Define Main System Parameters ---
P_PV_rated = 5000;
V_PV_mpp = 250;
I_PV_mpp = P_PV_rated/V_PV_mpp;
V_PV_oc = 300;
I_PV_sc = 22;
T_ref = 25;
G_ref = 1000;
P_Elec_rated = 4000;
P_fc_rated = 3000;
V_dc = 400;
% --- Step 5: Calculate Power Profile from Irradiance ---
effective_area_efficiency = P_PV_rated / G_ref;
solar_power_watts = solar_irradiance * effective_area_efficiency;
solar_power_kw = solar_power_watts / 1000;
solar_profile_power = [solar_profile_time, solar_power_kw];
% Load Power Profile
drive_cycle_power_kw = 2 * ones(size(drive_cycle_time)); % Creates a constant 2 kW load
drive_cycle_power = [drive_cycle_time, drive_cycle_power_kw];
% --- Step 6: Define Physical and Environmental Parameters ---
stack_num_cells = 50;
stack_area = 0.1;
stack_w_channels = 0.001;
stack_num_channels = 50;
water_pipe_D = 0.02;
gas_pipe_D = 0.02;
env_p = 101325 + 1e-6;    % Initial environmental pressure (Pa), with a small offset
env_T = 25;                 % Initial environmental temperature (degC)
env_RH = 0.5;               % Environmental relative humidity (0 to 1)
env_yO2 = 0.21;             % Mole fraction of O2 in the environment
cathode_tube_D = 0.02;					  
anode_tube_D = 0.02;
tank_yH2 = 0; 
tank_V = 0.1; 
tank_p = 1; 
stack_iL = 10;              % Stack current limit or reference (A)
% --- Step 7: Define Heat Exchanger and Radiator Parameters ---
exchanger_air_area_primary = 1;
exchanger_eta_fin = 0.9;
exchanger_air_area_fins = 5;
exchanger_L = 0.5;
exchanger_W = 0.2;
exchanger_tube_Leq = 0.1;
exchanger_tube_H = 0.01;
exchanger_N_tubes = 20;
% --- Step 8: Define Cooling System and Radiator Parameters ---
coolant_num_layers = 2;
coolant_num_passes = 4;
coolant_tube_D = 0.01;
coolant_w_channels = 0.002;
radiator_air_area_primary = 1;
radiator_air_area_fins = 5;
radiator_cp = 900;
radiator_eta_fin = 0.9;
radiator_L = 0.6;
radiator_W = 0.4;
radiator_tube_H = 0.01;
radiator_tube_Leq = 0.2;
radiator_N_tubes = 40;
radiator_t_wall = 0.001;
radiator_rho = 2700;
% --- Step 9: Define Compressor Map Parameters (Placeholders) ---
comp_rpm_TLU = [1000, 2000, 3000]; 
comp_p_ratio_TLU = [1.0, 1.5, 2.0];
comp_mdot_corr_TLU = [ 0.01, 0.009, 0.008; ...
    0.02, 0.018, 0.016; ...
    0.03, 0.027, 0.024 ];
disp('All parameters successfully loaded.');

Import:

Rodney Tan (2025). MPPT Solar Charge Controller Model (https://www.mathworks.com/matlabcentral/fileexchange/73115-mppt-solar-charge-controller-model), MATLAB Central File Exchange. Retrieved September 1, 2025.
PEM Electrolysis System > openExample('simscape/PEMElectrolysisSystemExample')

Link > PEM Electrolysis System

I've tried changing linear algebra setting in solver configuration to full/sparse/auto and using ode23t (mod. stiff/Trapezoidal) solver - Matlab 2024b. Any tips/ideas/clues?

Answer 1

[JAY]

Canvas layout - attached.