%% verify_deadtime_formula3 % Verify the dead-time voltage-error model behind formula (3) in % electronics-08-00196-v2.pdf. % % The script uses a normalized single inverter leg: % q = 0 means the phase node is clamped to the lower DC rail. % q = 1 means the phase node is clamped to the upper DC rail. % % Dead-time average error model: % delta_q = polarity * rho * sign(i) % % where rho = Td / Tc = fc * Td. clear; close all; clc; scriptFolder = fileparts(mfilename("fullpath")); baseCfg.name = "baseline_50Hz_10kHz"; baseCfg.description = "Baseline: 50 Hz fundamental, 10 kHz carrier, rho = 0.02"; baseCfg.fm = 50; % Fundamental/modulating frequency [Hz] baseCfg.fc = 10000; % Carrier frequency [Hz] baseCfg.M = 0.8; % Modulation index for d(t)=0.5+M/2*sin(wt) baseCfg.rho = 0.02; % Dead-time ratio, Td/Tc baseCfg.dead_time_s = baseCfg.rho / baseCfg.fc; baseCfg.power_factor = 0.8; % PF = cos(theta) baseCfg.current_lagging = true; % true: current lags the modulation voltage baseCfg.polarity = 1; % +1 or -1; selects the error sign convention baseCfg.switching_current_mode = "carrier-center"; baseCfg.fundamental_cycles = 8; % Coherent cycles for harmonic projection baseCfg.samples_per_fundamental = 20000; baseCfg.switch_cycles = 2; % Fundamental cycles in switching simulation baseCfg.samples_per_carrier = 400; baseCfg.max_harmonic = 15; baseCfg.figure_file = "deadtime_formula3_verification.png"; baseCfg.avg_waveform_file = "deadtime_formula3_avg_waveform_baseline.png"; baseCfg.harmonic_table_file = "deadtime_formula3_harmonics.csv"; aggressiveCfg = baseCfg; aggressiveCfg.name = "aggressive_1000Hz_20kHz_2us"; aggressiveCfg.description = "Aggressive: 1000 Hz fundamental, 20 kHz carrier, 2 us dead time"; aggressiveCfg.fm = 1000; aggressiveCfg.fc = 20000; aggressiveCfg.dead_time_s = 2e-6; aggressiveCfg.rho = aggressiveCfg.fc * aggressiveCfg.dead_time_s; aggressiveCfg.switching_current_mode = "instantaneous"; aggressiveCfg.switch_cycles = 4; aggressiveCfg.figure_file = "deadtime_formula3_verification_aggressive.png"; aggressiveCfg.avg_waveform_file = "deadtime_formula3_avg_waveform_aggressive.png"; aggressiveCfg.harmonic_table_file = "deadtime_formula3_harmonics_aggressive.csv"; scenarios = [baseCfg, aggressiveCfg]; summaries = repmat(emptySummary(), numel(scenarios), 1); fprintf("Dead-time formula (3) verification\n"); for idx = 1:numel(scenarios) summaries(idx) = runScenario(scenarios(idx), scriptFolder); end summaryTable = struct2table(summaries); summaryPath = fullfile(scriptFolder, "deadtime_formula3_summary.csv"); writetable(summaryTable, summaryPath); disp("Scenario summary:"); disp(summaryTable); fprintf("Saved scenario summary: %s\n", summaryPath); function summary = emptySummary() summary = struct( ... "scenario", "", ... "fm_Hz", NaN, ... "fc_Hz", NaN, ... "carrier_ratio", NaN, ... "dead_time_us", NaN, ... "rho", NaN, ... "power_factor", NaN, ... "theta_deg", NaN, ... "switching_current_mode", "", ... "fundamental_error_amplitude", NaN, ... "max_harmonic_amplitude_error", NaN, ... "max_phase_error_deg", NaN, ... "switching_rms_error", NaN, ... "switching_max_error", NaN); end function summary = runScenario(cfg, scriptFolder) validateConfig(cfg); theta = acos(cfg.power_factor); if ~cfg.current_lagging theta = -theta; end Tc = 1 / cfg.fc; Td = cfg.rho * Tc; carrierRatio = cfg.fc / cfg.fm; fprintf("\n=== %s ===\n", cfg.description); fprintf("fm = %.6g Hz, fc = %.6g Hz, fc/fm = %.3f\n", ... cfg.fm, cfg.fc, carrierRatio); fprintf("Tc = %.6g us, Td = %.6g us, rho = %.6g\n", ... Tc * 1e6, Td * 1e6, cfg.rho); fprintf("M = %.4f, power factor = %.4f, theta = %.4f deg, polarity = %+d\n", ... cfg.M, cfg.power_factor, rad2deg(theta), cfg.polarity); fprintf("Average-model settings: %d fundamental cycles, %d samples/cycle\n", ... cfg.fundamental_cycles, cfg.samples_per_fundamental); fprintf("Switching-model settings: %d fundamental cycles, %d samples/carrier, current mode = %s\n", ... cfg.switch_cycles, cfg.samples_per_carrier, cfg.switching_current_mode); fprintf("Max harmonic checked = %d\n\n", cfg.max_harmonic); [tAvg, modulation, current, deltaAvg] = buildAverageModel(cfg, theta); harmonics = oddHarmonicTable(tAvg, deltaAvg, cfg, theta); disp("Harmonic projection of delta_q = polarity*rho*sign(i):"); disp(harmonics); amplitudeError = harmonics.amplitude - harmonics.expected_amplitude; maxAmplitudeError = max(abs(amplitudeError)); maxPhaseErrorDeg = max(abs(harmonics.phase_error_deg)); fprintf("Max harmonic amplitude error = %.3e\n", maxAmplitudeError); fprintf("Max harmonic phase error = %.3e deg\n\n", maxPhaseErrorDeg); switching = simulateSwitchingDeadTime(cfg, theta, cfg.switching_current_mode); cycleError = switching.cycle_error - switching.expected_cycle_error; switchRmsError = sqrt(mean(cycleError.^2)); switchMaxError = max(abs(cycleError)); fprintf("Switching-model cycle-average RMS error = %.3e\n", switchRmsError); fprintf("Switching-model cycle-average max error = %.3e\n\n", switchMaxError); figurePath = fullfile(scriptFolder, char(cfg.figure_file)); avgWaveformPath = fullfile(scriptFolder, char(cfg.avg_waveform_file)); harmonicTablePath = fullfile(scriptFolder, char(cfg.harmonic_table_file)); writetable(harmonics, harmonicTablePath); plotVerificationFigure(tAvg, modulation, current, deltaAvg, switching, ... cfg, theta, figurePath); plotAverageOutputComparison(tAvg, modulation, current, deltaAvg, ... cfg, theta, avgWaveformPath); fprintf("Saved harmonic table: %s\n", harmonicTablePath); fprintf("Saved verification figure: %s\n", figurePath); fprintf("Saved average waveform comparison: %s\n", avgWaveformPath); summary = emptySummary(); summary.scenario = char(cfg.name); summary.fm_Hz = cfg.fm; summary.fc_Hz = cfg.fc; summary.carrier_ratio = carrierRatio; summary.dead_time_us = Td * 1e6; summary.rho = cfg.rho; summary.power_factor = cfg.power_factor; summary.theta_deg = rad2deg(theta); summary.switching_current_mode = char(cfg.switching_current_mode); summary.fundamental_error_amplitude = harmonics.expected_amplitude(1); summary.max_harmonic_amplitude_error = maxAmplitudeError; summary.max_phase_error_deg = maxPhaseErrorDeg; summary.switching_rms_error = switchRmsError; summary.switching_max_error = switchMaxError; end function validateConfig(cfg) arguments cfg struct end mustBePositiveScalar(cfg.fm, "cfg.fm"); mustBePositiveScalar(cfg.fc, "cfg.fc"); mustBePositiveScalar(cfg.dead_time_s, "cfg.dead_time_s"); mustBePositiveScalar(cfg.samples_per_fundamental, ... "cfg.samples_per_fundamental"); mustBePositiveScalar(cfg.samples_per_carrier, ... "cfg.samples_per_carrier"); if cfg.M < 0 || cfg.M > 1 error("verifyDeadtime:invalidM", ... "cfg.M must be between 0 and 1."); end if cfg.rho <= 0 || cfg.rho >= 0.5 error("verifyDeadtime:invalidRho", ... "cfg.rho must be in the open interval (0, 0.5)."); end if cfg.power_factor < 0 || cfg.power_factor > 1 error("verifyDeadtime:invalidPowerFactor", ... "cfg.power_factor must be between 0 and 1."); end if ~ismember(cfg.polarity, [-1, 1]) error("verifyDeadtime:invalidPolarity", ... "cfg.polarity must be +1 or -1."); end if ~ismember(string(cfg.switching_current_mode), ... ["carrier-center", "instantaneous"]) error("verifyDeadtime:invalidSwitchingCurrentMode", ... "cfg.switching_current_mode must be carrier-center or instantaneous."); end rhoFromDeadTime = cfg.fc * cfg.dead_time_s; if abs(rhoFromDeadTime - cfg.rho) > 1e-12 error("verifyDeadtime:inconsistentDeadTime", ... "cfg.rho must equal cfg.fc * cfg.dead_time_s."); end minPulseMargin = (1 - cfg.M) / 2; if cfg.rho >= minPulseMargin warning("verifyDeadtime:shortPulse", ... ["cfg.rho is not smaller than the minimum high/low pulse " ... "margin. The simple non-overlap switching model can break " ... "near the modulation extrema."]); end if cfg.fc <= 10 * cfg.fm warning("verifyDeadtime:lowCarrierRatio", ... "fc/fm is low; current polarity can change noticeably inside a carrier period."); end end function mustBePositiveScalar(value, name) if ~isscalar(value) || ~isnumeric(value) || value <= 0 error("verifyDeadtime:invalidPositiveScalar", ... "%s must be a positive numeric scalar.", name); end end function [t, modulation, current, deltaAvg] = buildAverageModel(cfg, theta) samples = cfg.fundamental_cycles * cfg.samples_per_fundamental; dt = 1 / (cfg.fm * cfg.samples_per_fundamental); t = (0:samples-1).' * dt; wt = 2 * pi * cfg.fm * t; modulation = 0.5 + 0.5 * cfg.M * sin(wt); current = sin(wt - theta); deltaAvg = cfg.polarity * cfg.rho * signNonzero(current); end function harmonics = oddHarmonicTable(t, signal, cfg, theta) harmonicNumber = (1:2:cfg.max_harmonic).'; amplitude = zeros(size(harmonicNumber)); phase_deg = zeros(size(harmonicNumber)); expected_amplitude = zeros(size(harmonicNumber)); expected_phase_deg = zeros(size(harmonicNumber)); signal = signal(:); t = t(:); for idx = 1:numel(harmonicNumber) n = harmonicNumber(idx); angle = 2 * pi * n * cfg.fm * t; cosCoeff = 2 / numel(t) * sum(signal .* cos(angle)); sinCoeff = 2 / numel(t) * sum(signal .* sin(angle)); amplitude(idx) = hypot(cosCoeff, sinCoeff); phase_deg(idx) = wrapTo180Local(atan2d(cosCoeff, sinCoeff)); expected_amplitude(idx) = 4 * cfg.rho / (pi * n); expected_phase_deg(idx) = wrapTo180Local(-n * rad2deg(theta)); if cfg.polarity < 0 expected_phase_deg(idx) = wrapTo180Local( ... expected_phase_deg(idx) + 180); end end amplitude_error = amplitude - expected_amplitude; phase_error_deg = wrapTo180Local(phase_deg - expected_phase_deg); harmonics = table(harmonicNumber, amplitude, expected_amplitude, ... amplitude_error, phase_deg, expected_phase_deg, phase_error_deg); end function switching = simulateSwitchingDeadTime(cfg, theta, signMode) Tc = 1 / cfg.fc; Td = cfg.rho * Tc; samplesPerCarrier = round(cfg.samples_per_carrier); carrierCount = max(1, round(cfg.switch_cycles * cfg.fc / cfg.fm)); totalSamples = carrierCount * samplesPerCarrier; dt = Tc / samplesPerCarrier; t = (0:totalSamples-1).' * dt; carrierIndex = floor(t / Tc) + 1; localTime = t - (carrierIndex - 1) * Tc; carrierCenterTime = ((0:carrierCount-1).' + 0.5) * Tc; duty = 0.5 + 0.5 * cfg.M * sin(2 * pi * cfg.fm * carrierCenterTime); currentAtCenter = sin(2 * pi * cfg.fm * carrierCenterTime - theta); currentSign = signNonzero(currentAtCenter); dutyAtSample = duty(carrierIndex); risingEdge = 0.5 * Tc - 0.5 * dutyAtSample * Tc; fallingEdge = 0.5 * Tc + 0.5 * dutyAtSample * Tc; idealHigh = localTime >= risingEdge & localTime < fallingEdge; inDeadTimeAfterRising = localTime >= risingEdge & localTime < risingEdge + Td; inDeadTimeAfterFalling = localTime >= fallingEdge & localTime < fallingEdge + Td; inDeadTime = inDeadTimeAfterRising | inDeadTimeAfterFalling; switch string(signMode) case "carrier-center" clampHigh = cfg.polarity * currentSign > 0; clampHighAtSample = clampHigh(carrierIndex); case "instantaneous" currentAtSample = sin(2 * pi * cfg.fm * t - theta); currentSignAtSample = signNonzero(currentAtSample); clampHighAtSample = cfg.polarity * currentSignAtSample > 0; otherwise error("verifyDeadtime:invalidSwitchingCurrentMode", ... "Unknown switching current mode."); end actualHigh = idealHigh; actualHigh(inDeadTime) = clampHighAtSample(inDeadTime); instantaneousError = double(actualHigh) - double(idealHigh); cycleError = mean(reshape(instantaneousError, samplesPerCarrier, []), 1).'; expectedCycleError = cfg.polarity * cfg.rho * currentSign; switching.t = t; switching.ideal_high = idealHigh; switching.actual_high = actualHigh; switching.instantaneous_error = instantaneousError; switching.carrier_center_time = carrierCenterTime; switching.duty = duty; switching.current_sign = currentSign; switching.current_mode = string(signMode); switching.cycle_error = cycleError; switching.expected_cycle_error = expectedCycleError; end function plotVerificationFigure(tAvg, modulation, current, deltaAvg, switching, ... cfg, theta, figurePath) oneCycle = tAvg < 1 / cfg.fm; reconstruction = reconstructDeadTimeError(tAvg, cfg, theta); fig = figure("Visible", "off", "Color", "w", ... "Position", [100, 100, 1050, 780]); subplot(3, 1, 1); yyaxis left; plot(tAvg(oneCycle) * 1e3, modulation(oneCycle), "LineWidth", 1.3); ylabel("m(t)"); ylim([0, 1]); yyaxis right; plot(tAvg(oneCycle) * 1e3, current(oneCycle), "LineWidth", 1.0); ylabel("i(t), normalized"); grid on; title(sprintf("%s, PF = %.3f, theta = %.2f deg", ... cfg.name, cfg.power_factor, rad2deg(theta)), "Interpreter", "none"); subplot(3, 1, 2); plot(tAvg(oneCycle) * 1e3, deltaAvg(oneCycle), "k", "LineWidth", 1.2); hold on; plot(tAvg(oneCycle) * 1e3, reconstruction(oneCycle), "r--", ... "LineWidth", 1.0); grid on; ylabel("\Delta q"); legend("sign(i) average model", ... sprintf("Fourier reconstruction to %dth", cfg.max_harmonic), ... "Location", "best"); subplot(3, 1, 3); plot(switching.carrier_center_time * 1e3, switching.cycle_error, ... "bo", "MarkerSize", 3); hold on; plot(switching.carrier_center_time * 1e3, ... switching.expected_cycle_error, "r.", "MarkerSize", 8); grid on; xlabel("time [ms]"); ylabel("cycle-avg \Delta q"); legend(sprintf("measured switching waveform (%s)", switching.current_mode), ... "polarity*rho*sign(i) at carrier center", ... "Location", "best"); exportgraphics(fig, figurePath, "Resolution", 160); close(fig); end function plotAverageOutputComparison(t, modulation, current, deltaAvg, cfg, ... theta, figurePath) oneCycle = t < 1 / cfg.fm; tMs = t(oneCycle) * 1e3; modulationOneCycle = modulation(oneCycle); actualAverageOneCycle = modulation(oneCycle) + deltaAvg(oneCycle); deltaOneCycle = deltaAvg(oneCycle); currentOneCycle = current(oneCycle); zeroCrossingsMs = currentZeroCrossings(t(oneCycle), currentOneCycle) * 1e3; fig = figure("Visible", "off", "Color", "w", ... "Position", [100, 100, 1050, 720]); subplot(2, 1, 1); plot(tMs, modulationOneCycle, "b", "LineWidth", 1.4); hold on; plot(tMs, actualAverageOneCycle, "r", "LineWidth", 1.4); drawZeroCrossingLines(zeroCrossingsMs); grid on; ylim([0, 1]); ylabel("cycle-average q"); title(sprintf("%s: ideal m(t) versus m(t)+\\Delta q, rho = %.3f, theta = %.2f deg", ... cfg.name, cfg.rho, rad2deg(theta)), "Interpreter", "none"); legend("ideal m(t)", "with dead-time average m(t)+\Delta q", ... "current zero crossing", "Location", "best"); subplot(2, 1, 2); yyaxis left; stairs(tMs, deltaOneCycle, "k", "LineWidth", 1.3); ylabel("\Delta q"); ylim(1.35 * cfg.rho * [-1, 1]); yyaxis right; plot(tMs, currentOneCycle, "Color", [0.85, 0.33, 0.1], "LineWidth", 1.0); ylabel("i(t), normalized"); ylim([-1.1, 1.1]); drawZeroCrossingLines(zeroCrossingsMs); grid on; xlabel("time [ms]"); title("Dead-time average error flips with current polarity"); legend("\Delta q = \rho sign(i)", "current", ... "current zero crossing", "Location", "best"); exportgraphics(fig, figurePath, "Resolution", 160); close(fig); end function zeroCrossings = currentZeroCrossings(t, current) signCurrent = signNonzero(current); crossingIndex = find(diff(signCurrent) ~= 0); zeroCrossings = zeros(size(crossingIndex)); for idx = 1:numel(crossingIndex) k = crossingIndex(idx); t1 = t(k); t2 = t(k + 1); i1 = current(k); i2 = current(k + 1); zeroCrossings(idx) = t1 - i1 * (t2 - t1) / (i2 - i1); end end function drawZeroCrossingLines(zeroCrossingsMs) for idx = 1:numel(zeroCrossingsMs) if idx == 1 xline(zeroCrossingsMs(idx), "k:", "current zero crossing", ... "LabelOrientation", "horizontal"); else xline(zeroCrossingsMs(idx), "k:", "HandleVisibility", "off"); end end end function y = reconstructDeadTimeError(t, cfg, theta) y = zeros(size(t)); for n = 1:2:cfg.max_harmonic y = y + cfg.polarity * 4 * cfg.rho / (pi * n) .* ... sin(n * (2 * pi * cfg.fm * t - theta)); end end function s = signNonzero(x) s = ones(size(x)); s(x < 0) = -1; end function angle = wrapTo180Local(angle) angle = mod(angle + 180, 360) - 180; end