Skip to content

Conversation

@simmsa
Copy link
Contributor

@simmsa simmsa commented Mar 12, 2025

This PR removes the python calls to MHKiT-Python and rewrites the MHKiT-Python functionality in native MATLAB code:

  • mhkit/power/characteristics/ac_three_phase.m
  • mhkit/power/characteristics/dc_power.m
  • mhkit/power/characteristics/instantaneous_frequency.m
  • mhkit/power/quality/calc_electrical_angle.m (Already native matlab, unchanged)
  • mhkit/power/quality/calc_flicker_coefficient.m (Already native matlab, unchanged)
  • mhkit/power/quality/calc_fundamental_freq.m (Already native matlab, unchanged)
  • mhkit/power/quality/calc_ideal_voltage.m (Already native matlab, unchanged)
  • mhkit/power/quality/calc_Rfic_Lfic.m (Already native matlab, unchanged)
  • mhkit/power/quality/calc_shortterm_flicker_severity.m (Already native matlab, unchanged)
  • mhkit/power/quality/calc_simulated_voltage.m (Already native matlab, unchanged)
  • mhkit/power/quality/flicker_ufic_workflow.m (Already native matlab, unchanged)
  • mhkit/power/quality/gen_test_data.m (Already native matlab, unchanged)
  • mhkit/power/quality/harmonic_subgroups.m
  • mhkit/power/quality/harmonics.m
  • mhkit/power/quality/interharmonics.m
  • mhkit/power/quality/total_harmonic_current_distortion.m

@simmsa simmsa changed the title Power Module: Convert module to native MATLAB code Power Module: Convert to native MATLAB code Mar 12, 2025
@simmsa simmsa requested review from hivanov-nrel and rpauly18 July 1, 2025 15:43
@simmsa
Copy link
Contributor Author

simmsa commented Jul 1, 2025

@rpauly18 and @hivanov-nrel , this PR is ready for review. I'm keeping this as a draft for now because I may need to make other updates to the tests for some other PR's and merge those into this PR.

Happy to answer any questions.

@simmsa simmsa requested a review from akeeste July 7, 2025 15:34
Copy link
Contributor

@akeeste akeeste left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @simmsa. This module is not a strength of mine, but everything seems to line up with the python implementation. The code is all very clean, well documented, and the function argument validation looks robust

end

% Validate time vectors match (if both are structures)
if isstruct(voltage) && isstruct(current)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't the time vectors match regardless of whether they come from a structure or matrix?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, good catch. This check can be removed.

if grid_freq == 60
max_freq = 3060; % Exactly 51st harmonic (60 * 51)
elseif grid_freq == 50
max_freq = 2570; % 51.4th harmonic (intentional extension)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why 51.4?

Copy link
Contributor Author

@simmsa simmsa Jul 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question, I am not entirely sure. This is from the python implementation: https://github.com/MHKiT-Software/MHKiT-Python/blob/68fe393855d4f8c71f2c6f15feada8464a2bc99e/mhkit/power/quality.py#L109. I need to look into the 61000-4-7 standard.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rpauly18, do you have any insight here. The intent of the harmonics code seems to be to calculate the input values to downstream harmonics quantities of interest calculations including total harmonic distortion:

image

Ref: IEC 61000-4-7 Edition 2.1 2009-10 Equation 4

In the code section referenced here $h_{max}$ is 50 and max_freq should be 2550. Does this seem correct to you?

Should we also add an argument so the user can set $h_{max}$ and verify the value in within some reasonable range, maybe 40 to 100?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am looking into this, but I am having a hard time following the history of this function on the Python side to figure out why this was developed this way. This function was originally different. The 51 Hz thing is very weird. The commit history on this function is also weird and inconsistent.

Comment on lines +78 to +80
if any(abs(time_diff - dt_mean) > dt_tolerance)
warning('MHKiT:instantaneous_frequency: Non-uniform time spacing detected, using local time differences');
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this warning necessary? The function will always use the local time difference rather than the dt_mean

% Calculate power for each phase
if line_to_line
% For line-to-line measurements, apply sqrt(3) correction
power_per_phase = abs_current .* (abs_voltage * sqrt(3));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

alternatively, you could send the data to dc_power(), take the gross power output and multiply the power factor and line_to_line factor

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is a good suggestion, but might be confusing to a user (Why are we using the dc power to calculate ac_power), even though the equations are the same. The python implementation does not call dc_power: https://github.com/MHKiT-Software/MHKiT-Python/blob/68fe393855d4f8c71f2c6f15feada8464a2bc99e/mhkit/power/characteristics.py#L247, so maybe it makes the most sense to use the same pattern here.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@simmsa The python code calls dc_power() a few lines above: https://github.com/MHKiT-Software/MHKiT-Python/blob/68fe393855d4f8c71f2c6f15feada8464a2bc99e/mhkit/power/characteristics.py#L239. The current implementation is fine as is though

simmsa added 2 commits July 17, 2025 14:49
This is redundant based on the checks performed before this check.
@simmsa simmsa added this to the 1.1 milestone Sep 30, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants