预览截图
应用介绍
免费下载Matlab 电能质量仿真工具，带有GUI界面
%% Power Quality Tool
% A lightweight application for investigating common power quality measures of
% arbitrary power system measurement waveforms and groups of such waveforms.
function PowerQualityTool
% Create invisible figure window with no title, menus, or toolbar.
f = figure('Visible', 'off', ...
'NumberTitle', 'off', ...
'MenuBar', 'none', ...
'ToolBar', 'none', ...
'Position', [360, 500, 460, 315]);
% Create tabs for time and frequency domain plots for easy switching.
hPlotTabGroup = uitabgroup('Parent', f, ...
'Units', 'Pixels', ...
'Position', [0, 0, 305, 315], ...
'SelectionChangedFcn', {@tabSelectCallback});
hTimePlotTab = uitab('Parent', hPlotTabGroup, ...
'Title', 'Time Domain');
hFrequencyPlotTab = uitab('Parent', hPlotTabGroup, ...
'Title', 'Frequency Domain');
% Create buttons, label, and drop-down, and then align them.
hLoad = uicontrol('Style', 'pushbutton', ...
'String', 'Load', ...
'Position', [345, 250, 70, 25], ...
'Callback', {@loadButtonCallback});
hLabel = uicontrol('Style', 'text', ...
'String', 'Select Data', ...
'Position', [355, 225, 60, 15]);
hPopup = uicontrol('Style', 'popupmenu', ...
'String', {'No Data'}, ...
'Position', [330, 195, 100, 25], ...
'Callback', {@popupMenuCallback});
hPanel = uipanel('Parent', f, ...
'Title', 'Results', ...
'Units', 'Pixels', ...
'Position', [315, 10, 125, 170]);
hTable = uitable('Parent', hPanel, ...
'Data', {'<HTML>f<sub>0</sub></HTML>', 0; ...
'<HTML>V<sub>RMS</sub></HTML>', 0; ...
'<HTML>I<sub>RMS</sub></HTML>', 0; ...
'<HTML>THD<sub>V</sub></HTML>', 0; ...
'<HTML>THD<sub>I</sub></HTML>', 0; ...
'PF', 0; ...
'DF', 0}, ...
'RowName', {}, ...
'ColumnName', {'', 'Value'}, ...
'ColumnWidth', {35, 'auto'}, ...
'Position', [5, 5, 115, 150]);
align([hLoad, hLabel, hPopup, hPanel, hTable], 'Center', 'None');
% Create axes for both time and frequency plots within respective tabs.
hTimeAxes = axes('Parent', hTimePlotTab, ...
'Units', 'Pixels', ...
'Position', [50, 60, 200, 200]);
hFrequencyAxes = axes('Parent', hFrequencyPlotTab, ...
'Units', 'Pixels', ...
'Position', [50, 60, 200, 200]);
% Change to normalized units for relative positioning (vs. absolute).
f.Units = 'normalized';
hPlotTabGroup.Units = 'normalized';
hTimePlotTab.Units = 'normalized';
hFrequencyPlotTab.Units = 'normalized';
hTimeAxes.Units = 'normalized';
hFrequencyAxes.Units = 'normalized';
hLoad.Units = 'normalized';
hLabel.Units = 'normalized';
hPopup.Units = 'normalized';
hPanel.Units = 'normalized';
hTable.Units = 'normalized';
% Assign a name to the tool's title bar.
f.Name = 'Power Quality Tool';
% Center the window on the screen.
movegui(f, 'center')
% Show window with axes and controls.
f.Visible = 'on';
% Create data and results structures.
Data = struct(); % for user input data
Results = struct(); % for calculated values
currentData = [];
currentName = '';
currentTab = hPlotTabGroup.SelectedTab.Title;
% Callback for the tab group tab selection.
function tabSelectCallback(source, eventdata)
currentTab = eventdata.NewValue.Title;
if ~strcmp(currentName, '')
updatePlot(currentData, currentName, currentTab, Results);
end
end
% Callback for the drop-down menu.
function popupMenuCallback(source, eventdata)
% Determine if single or multiple data sets (buses) and assign a set.
[~, numData] = size(Data);
currentName = source.String{source.Value};
if numData == 1
currentData = [Data.Time, Data.Voltage, Data.Current];
else
index = find(strcmp(currentName, {Data.Name}));
currentData = [Data(index).Time, Data(index).Voltage, Data(index).Current];
end
% Update results data and plot of selected data set.
[currentData, Results] = updateResults(currentData, Results);
updatePlot(currentData, currentName, currentTab, Results);
end
% Callback for 'Load' button. Requires data in CSV format, with headers
% and arranged in columns: Time, Voltage, Current.
function loadButtonCallback(source, eventdata)
% Let the user pick a (few) file(s).
fileTypesCell = {'*.csv;*.txt;*.dat', ...
'Comma-Separated Values (*.csv,*.txt,*.dat)'; ...
'*.*', ...
'All Files (*.*)'};
[fileName, pathName] = uigetfile(fileTypesCell, ...
'Select waveform data', ...
'MultiSelect', 'on');
% Read and prepare the data.
if ~(isequal(fileName, 0) || isequal(pathName, 0)) % user didn't cancel
if iscell(fileName)
[~, cellSize] = size(fileName);
for i = 1:cellSize
tempData = csvread([pathName, fileName{i}], 1, 0);
Data(i).Name = fileName{i}(1:end - 4);
Data(i).Time = tempData(:, 1) - tempData(1, 1); % remove offset
Data(i).Voltage = tempData(:, 2);
Data(i).Current = tempData(:, 3);
end
currentData = [Data(1).Time, Data(1).Voltage, Data(1).Current];
else
tempData = csvread([pathName, fileName], 1, 0);
Data.Name = fileName(1:end - 4);
Data.Time = tempData(:, 1) - tempData(1, 1); % remove offset
Data.Voltage = tempData(:, 2);
Data.Current = tempData(:, 3);
currentData = [Data.Time, Data.Voltage, Data.Current];
end
% Inject the data and its identifiers into the appropriate places.
hPopup.String = {Data.Name};
% Automatically select the first file.
hPopup.Value = 1;
currentName = hPopup.String{1};
% Update result data and plot with (first) loaded data.
[currentData, Results] = updateResults(currentData, Results);
updatePlot(currentData, currentName, currentTab, Results);
end
end
% Result update function. Updates RMS, THD, PF and DF. Updates frequency
% domain values for both voltage and current.
function [currentData, Results] = updateResults(currentData, Results)
% Gather time information.
period60 = 1/60; % static 60 Hz cycle definition
Results.Tm = currentData(end, 1); % measurement period
% Create vertical line t parameters to indicate periods of 60 Hz.
Results.timeVerts = struct();
Results.nPeriods = floor(Results.Tm/period60);
for iVerts = 1:Results.nPeriods
Results.timeVerts(iVerts).t = [iVerts*period60, iVerts*period60];
end
[currentData, Results] = cleanData(currentData, Results);
% Calculate RMS voltage and current of composite waveform.
Results.RMSv = rms(currentData(:, 2));
Results.RMSi = rms(currentData(:, 3));
% Calculate and store peaks for time domain annotations.
[Results.vtPeaks, Results.vtPeakIndices] = findpeaks(currentData(:, 2));
[Results.itPeaks, Results.itPeakIndices] = findpeaks(currentData(:, 3));
% Find fundamental frequency allowing search of entirety of waveforms.
Results = findFundamentalFrequency(currentData, Results);
fftVoltageResult = abs(fft(Results.correctedVoltage)/(Results.nTimePoints + 1));
fftVoltageResult((Results.nTimePoints - 1)/2 + 1:end) = [];
fftVoltageResult(2:end - 1) = 2*fftVoltageResult(2:end - 1);
fftCurrentResult = abs(fft(Results.correctedCurrent)/(Results.nTimePoints + 1));
fftCurrentResult((Results.nTimePoints - 1)/2 + 1:end) = [];
fftCurrentResult(2:end - 1) = 2*fftCurrentResult(2:end - 1);
% Truncate all results to 5 kHz.
Results.frequency = Results.sampleFrequency*(0:(Results.nTimePoints/2 + 1))/Results.nTimePoints;
Results.frequency(Results.frequency > 5000) = [];
Results.frequency = Results.frequency';
numFrequencies = length(Results.frequency);
Results.fVoltage = fftVoltageResult(1:numFrequencies);
Results.fCurrent = fftCurrentResult(1:numFrequencies);
% Calculate power factor.
dt = 1/Results.sampleFrequency;
Pavg = Results.correctedVoltage'*Results.correctedCurrent*dt/(Results.nPeriods/60);
Results.powerFactor = Pavg/(Results.RMSv*Results.RMSi);
% Calculate harmonic component amplitudes.
harmonicFrequencies = 60*(1:floor(5e3/60)); % harmonics to 5 kHz
Results.harmonicVoltage = [];
Results.harmonicCurrent = [];
for harmonic = 1:length(harmonicFrequencies)
[frequencyIndex, ~] = find(Results.frequency == harmonicFrequencies(harmonic));
Results.harmonicVoltage = [Results.harmonicVoltage; Results.fVoltage(frequencyIndex)];
Results.harmonicCurrent = [Results.harmonicCurrent; Results.fCurrent(frequencyIndex)];
end
% Calculate THDv and THDi.
Results.THDv = norm(Results.harmonicVoltage(2:end), 2)/Results.harmonicVoltage(1);
Results.THDi = norm(Results.harmonicCurrent(2:end), 2)/Results.harmonicCurrent(1);
% Calculate displacement factor using PF and THDi/v.
Results.displacementFactor = sqrt(1 + Results.THDv^2)*sqrt(1 + Results.THDi^2)*Results.powerFactor;
% Calculate and store peaks for frequency domain annotations.
prominenceV = 0.01*max(Results.fVoltage);
prominenceI = 0.01*max(Results.fCurrent);
[Results.vfPeaks, Results.vfPeakIndices] = findpeaks(Results.fVoltage, ...
'MinPeakProminence', prominenceV);
[Results.ifPeaks, Results.ifPeakIndices] = findpeaks(Results.fCurrent, ...
'MinPeakProminence', prominenceI);
% Assign data to table for display.
tableData = tablePrep(Results);
hTable.Data(:, 2) = tableData;
end
% Cleans data for processing by smoothing and extending length to a power of
% two number of data points via spline interpolation.
function [cleanedData, Results] = cleanData(currentData, Results)
% Define final measurement time point.
Tm = currentData(end, 1);
% Pick out last periodic point and set aside periodic segment.
sampleTimeGoal = Results.nPeriods/60;
[~, timeIndex] = min(abs(currentData(:, 1) - sampleTimeGoal));
Tp = currentData(timeIndex, 1);
truncationIndices = 1:timeIndex - 1;
timeVector = currentData(truncationIndices, 1);
nTimePoints = length(timeVector);
% Interpolate
nFullTimePoints = ceil(2^(nextpow2(nTimePoints))*Tm/Tp);
cleanedData(:, 1) = linspace(0, Tm, nFullTimePoints);
cleanedData(:, 2) = spline(currentData(:, 1), currentData(:, 2), cleanedData(:, 1));
cleanedData(:, 3) = spline(currentData(:, 1), currentData(:, 3), cleanedData(:, 1));
sampledVoltage = currentData(truncationIndices, 2);
sampledCurrent = currentData(truncationIndices, 3);
% Prepare for FFT.
nTimePointsGoal = 2^nextpow2(nTimePoints);
timeVectorGoal = linspace(0, sampleTimeGoal, nTimePointsGoal);
correctedVoltage = spline(timeVector, sampledVoltage, timeVectorGoal)';
correctedCurrent = spline(timeVector, sampledCurrent, timeVectorGoal)';
% Assign values for next segment
Results.timeVector = timeVectorGoal(1:end - 1);
Results.nTimePoints = length(Results.timeVector);
Results.sampledVoltage = sampledVoltage;
Results.correctedVoltage = correctedVoltage(1:end - 1);
Results.sampledCurrent = sampledCurrent;
Results.correctedCurrent = correctedCurrent(1:end - 1);
Results.sampleFrequency = Results.nTimePoints/timeVectorGoal(end);
end
% Fundamental frequency detection function. Finds fundamental frequency by
% averaging all distances to peaks within 0.1% magnitude of each other.
function Results = findFundamentalFrequency(currentData, Results)
time = currentData(:, 1);
dt = mean(diff(time));
freq = [];
vpks = Results.vtPeaks;
vpkt = Results.vtPeakIndices;
ipks = Results.itPeaks;
ipkt = Results.itPeakIndices;
% Periodic testing parameters.
runningAverage = 1/60; % initialize to line frequency
identified = 0;
% Find time gaps between peaks of similar voltage magnitude.
for vPeak = 1:length(vpks)
tempIdx = vpkt;
tempIdx(vPeak) = []; % exclude test point
vField = vpks;
vField(vPeak) = [];
vDiff = vpks(vPeak) - vField; % field points
for idx = 1:length(vDiff)
if vDiff(idx) <= 0.0001*vpks(vPeak) % difference is <= 0.01%
timeDiff = abs(time(tempIdx(idx)) - time(vpkt(vPeak)));
count1 = 1;
count2 = 1;
% Test for multiple/fractional period.
while ~identified
ratio = timeDiff/runningAverage;
if ratio > 1.01 % more than one period +/- 1%
count2 = count2 + 1;
timeDiff = timeDiff*count1/count2;
elseif ratio < 0.99 % less than one period +/- 1%
count1 = count1 + 1;
timeDiff = timeDiff*count1/count2;
else % about one period +/- 10%
if rem(count1, count2) == 0 % near integer period, keep
freq = [freq; 1/timeDiff];
runningAverage = mean([1./freq; runningAverage]);
identified = 1;
else % we have a fraction of a period, discard
identified = 1;
end
end
end
end
identified = 0;% prepare for next iteration
end
end
% Find time gaps between peaks of similar current magnitude.
for iPeak = 1:length(ipks)
tempIdx = ipkt;
tempIdx(iPeak) = []; % exclude test point
iField = ipks;
iField(iPeak) = [];
iDiff = ipks(iPeak) - iField; % field points
for idx = 1:length(iDiff)
if iDiff(idx) <= 0.0001*ipks(iPeak) % difference is <= 0.01%
timeDiff = abs(time(tempIdx(idx)) - time(ipkt(iPeak)));
count1 = 1;
count2 = 1;
% Test for multiple/fractional period.
while ~identified
ratio = timeDiff/runningAverage;
if ratio > 1.01 % more than one period + 1%
count2 = count2 + 1;
timeDiff = timeDiff*count1/count2;
elseif ratio < 0.99 % less than one period - 1%
count1 = count1 + 1;
timeDiff = timeDiff*count1/count2;
else % about one period +/- 1%
if rem(count1, count2) == 0 % near integer period, keep
freq = [freq; 1/timeDiff];
runningAverage = mean([1./freq; runningAverage]);
identified = 1;
else % we have a fraction of a period, discard
identified = 1;
end
end
end
identified = 0;% prepare for next iteration
end
end
end
Results.f0 = mean(freq);
end
% Table data preparation function. Formats data according to type,
% magnitude, and precision.
function tableData = tablePrep(Results)
tableData = {Results.f0; ...
Results.RMSv; ...
Results.RMSi; ...
Results.THDv; ...
Results.THDi; ...
Results.powerFactor; ...
Results.displacementFactor};
% Format fundamental frequency for precision.
tableData{1} = sprintf('%3.1f Hz', tableData{1});
% Format RMS voltage for expected ranges.
if Results.RMSv >= 1e6
tableData{2} = sprintf('%5.3f MV', tableData{2}/1e6);
elseif Results.RMSv >= 1e3
tableData{2} = sprintf('%5.3f kV', tableData{2}/1e3);
else
tableData{2} = sprintf('%5.3f V', tableData{2});
end
% Format RMS current for expected ranges.
if Results.RMSi >= 1e3
tableData{3} = sprintf('%5.3f kA', tableData{3}/1e3);
elseif Results.RMSi >= 1
tableData{3} = sprintf('%5.3f A', tableData{3});
else
tableData{3} = sprintf('%5.3f mA', tableData{3}/1e-3);
end
% Format voltage THD as a percent.
tableData{4} = sprintf('%4.4f %%', tableData{4}*100);
% Format current THD as a percent.
tableData{5} = sprintf('%4.4f %%', tableData{5}*100);
% Format power factor, keeping up to 4 decimal places.
tableData{6} = sprintf('%.4f', tableData{6});
% Format displacement factor, keeping only two decimal places.
tableData{7} = sprintf('%.2f', tableData{7});
end
% Plot update function. Updates appropriate plot with fresh data after load
% or on drop-down selection.
function updatePlot(currentData, currentName, selectedTab, Results)
switch selectedTab
case 'Time Domain'
axes(hTimeAxes);
cla reset;
yyaxis left
plot(currentData(:, 1), currentData(:, 2));
title([currentName, ' Waveforms'])
xlabel('Time')
ylabel('Voltage')
grid on;
yyaxis right
plot(currentData(:, 1), currentData(:, 3));
ylabel('Current')
% Include period indicator lines.
for iVerts = 1:Results.nPeriods
Results.timeVerts(iVerts).y = hTimeAxes.YLim;
x = Results.timeVerts(iVerts).t;
y = Results.timeVerts(iVerts).y;
line(x, y, 'Color', 'black');
end
% Annotate peaks.
ivVals = [max(Results.vtPeaks), max(Results.itPeaks)];
tVals = [currentData(Results.vtPeakIndices(Results.vtPeaks == ivVals(1)), 1), ...
currentData(Results.itPeakIndices(Results.itPeaks == ivVals(2)), 1)];
labels = {sprintf('V_{pk} = %.2g V', ivVals(1)), ...
sprintf('I_{pk} = %.2g A', ivVals(2))};
% Assign text location based on peak time location.
if tVals(1) <= currentData(end, 1)/2
directionV = 'left';
else
directionV = 'right';
end
if tVals(2) <= currentData(end, 1)/2
directionI = 'left';
else
directionI = 'right';
end
% Place voltage peak annotation.
yyaxis left
vPeakLabel = text(tVals(1), ivVals(1), labels{1});
vPeakLabel.VerticalAlignment = 'top';
vPeakLabel.HorizontalAlignment = directionV;
vPeakLabel.Color = hTimeAxes.YAxis(1).Color;
% Place current peak annotation.
yyaxis right
iPeakLabel = text(tVals(2), ivVals(2), labels{2});
iPeakLabel.VerticalAlignment = 'bottom';
iPeakLabel.HorizontalAlignment = directionI;
iPeakLabel.Color = hTimeAxes.YAxis(2).Color;
case 'Frequency Domain'
axes(hFrequencyAxes);
cla reset;
yyaxis left
bar(Results.frequency/60, Results.fVoltage, 6)
title([currentName, ' Spectrum'])
xlabel('Harmonic')
ylabel('Voltage')
grid on
yyaxis right
bar(Results.frequency/60, Results.fCurrent, 3)
ylabel('Current')
% Create significant (1%) harmonic annotations.
vVals = Results.vfPeaks;
iVals = Results.ifPeaks;
fValsV = Results.frequency(Results.vfPeakIndices)/60;
fValsI = Results.frequency(Results.ifPeakIndices)/60;
vLabels = [];
iLabels = [];
for peak = 1:length(vVals)
vLabels = [vLabels, sprintf('V_{%d} = %3.4g V\n', fValsV(peak), vVals(peak))];
end
for peak = 1:length(iVals)
iLabels = [iLabels, sprintf('I_{%d} = %3.4g A\n', fValsI(peak), iVals(peak))];
end
% Place peak annotations in legend.
legend(vLabels, iLabels)
otherwise
end
end
end
©版权声明：本文内容由互联网用户自发贡献，版权归原创作者所有，本站不拥有所有权，也不承担相关法律责任。如果您发现本站中有涉嫌抄袭的内容，欢迎发送邮件至： [email protected] 进行举报，并提供相关证据，一经查实，本站将立刻删除涉嫌侵权内容。
转载请注明出处： apollocode » 电能质量工具
文件列表(部分)
| 名称 | 大小 | 修改日期 |
|---|---|---|
| license.txt | 0.70 KB | 2017-11-21 |
| PowerQualityTool.m | 5.10 KB | 2017-11-21 |
| PowerQuality_powermatlab | 0.00 KB | 2018-03-22 |
发表评论 取消回复