电力系统最优潮流问题（ＯＰＦ）可描述为在满足电力系统各种运行约束限制前提条件下，通过优化得到电力网络中各控制变量取值，使系统的某个性能指标达到最优，因此ＯＰＦ问题是一个典型的约束非线性规划问题。以简化梯度法、牛顿法、内点法、解耦法为代表的传统优化方法已被运用于解决ＯＰＦ问题上，但简化梯度法会出现锯齿现象，收敛性较差，且计算量大耗时多；要求目标函数具有光滑、连续可微的特性限制了牛顿法在解决ＯＰＦ问题中的应用；内点法计算速度慢、易收敛至局部极值。针对经典优化算法存在的不足，帝企鹅启发式智能优化算法，为解决ＯＰＦ问题的一种新途径［１］
function om = opf_model(mpc)
%OPF_MODEL Constructor for OPF model class.
% OM = OPF_MODEL(MPC)
%
% This class implements the OPF model object used to encapsulate
% a given OPF problem formulation. It allows for access to optimization
% variables, constraints and costs in named blocks, keeping track of the
% ordering and indexing of the blocks as variables, constraints and costs
% are added to the problem.
%
% This class is a sub-class of OPT_MODEL and simply adds the 'mpc'
% field for storing the MATPOWER case struct used to build the object
% along with the get_mpc() method.
%
% The following is the structure of the data in the OPF model object.
% Each field of .idx or .data is a struct whose field names are the names
% of the corresponding blocks of vars, constraints or costs (found in
% order in the corresponding .order field). The description next to these
% fields gives the meaning of the value for each named sub-field.
% E.g. om.var.data.v0.Pg contains a vector of initial values for the 'Pg'
% block of variables.
%
% om
% .opt_model - the corresponding OPT_MODEL object
% .mpc - MATPOWER case struct used to create this model object
% .baseMVA
% .bus
% .branch
% .gen
% .gencost
% .A (if present, must have l, u)
% .l
% .u
% .N (if present, must have fparm, H, Cw)
% .fparm
% .H
% .Cw
%
% See also OPT_MODEL.
es = struct();
s = struct('mpc', es);
om = opt_model;
om = class(s, 'opf_model', om);
else
if isa(mpc,'opf_model')
om = mpc;
else
if isfield(mpc, 'om') %% avoid nesting
s = struct('mpc', rmfield(mpc, 'om'));
else
s = struct('mpc', mpc);
end
om = opt_model;
om = class(s, 'opf_model', om);
end
end
function results = maxloadlim(mpc,dir_mll,varargin)
% MAXLOADLIM computes the maximum loadability limit in one direction. It
% uses dispatchable loads in MATPOWER
% RESULTS = MAXLOADLIM(MPC,DIR_MLL) returns the results from the
% optimization problem looking for the maximum loadability limit in
% the direction of load increase DIR_MLL. DIR_MLL defines the directions
% of load increases for all buses. For buses with zero loads, the
% direction of load increases must be zero. RESULTS contains all fields
% returned from the runopf MATPOWER function. It also contains the
% following additional fields:
% * dir_mll: the direction of load increase used as input.
% * stab_marg: the stability margin to the maximum loadability point from
% the base case defined in the input MPC.
% * bif: information about the bifurcation at the MLL point.
%
% RESULTS = MAXLOADLIM(MPC,DIR_MLL,NAME,VALUE) uses the options defined
% by the pair NAME,VALUE. The currently supported options are
% * 'verbose': 1 or 0 (Default). If set to 1, a summary of the results
% at the maximum loadability limit is printed.
% * 'use_qlim': 1 (Default) or 0. Enforces or not the reactive power
% limits of the generators.
% * 'Vlims_bus_nb': [] (Default) or array of integers. By default, the
% bus voltage limits are not enforced. This option allows for defining
% a set of buses at which the voltage limits are enforced.
%
% See also PREPARE_MAXLOADLIM, POSTPROC_MAXLOADLIM, PRINT_MAXLOADLIM,
% RUNOPF.
% MATPOWER
% Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
% by Camille Hamon
%
% This file is part of MATPOWER.
% Covered by the 3-clause BSD License (see LICENSE file for details).
% See http://www.pserc.cornell.edu/matpower/ for more info.
define_constants;
%% Checking the options, if any
input_checker = inputParser;
default_verbose = 0;
verbose_levels = [0;1];
check_verbose = @(x)(isnumeric(x) && isscalar(x) && any(x == verbose_levels));
addParameter(input_checker,'verbose',default_verbose,check_verbose);
input_checker.KeepUnmatched = true;
parse(input_checker,varargin{:});
options = input_checker.Results;
%% Prepare the matpower case for the maximum loadability limit problem
mpc_vl = prepare_maxloadlim(mpc,dir_mll,varargin{:});
%% Run opf
% Turning off the printing and initializing from the base case
mpopt = mpoption('verbose',options.verbose,'opf.init_from_mpc',1);
mpopt = mpoption(mpopt,'out.all',0);
% Decreasing the threshold for the relative complementarity constraints
mpopt = mpoption(mpopt,'mips.comptol',1e-8);
% Change solver
mpopt = mpoption(mpopt, 'opf.ac.solver', 'MIPS');
% Execute opf
results = runopf(mpc_vl,mpopt);
%% Post-processing
results = postproc_maxloadlim(results,dir_mll);
%% Printing
if options.verbose
print_maxloadlim(mpc,results);
end
[1]李英. 基于并行计算和粒子群优化算法的电力系统无功优化问题研究[D]. 浙江大学, 2010.
[2]张东寅, 王澎涛, 袁艳斌,等. 基于改进布谷鸟算法的电力系统最优潮流计算[J]. 水电能源科学, 2017, 35(1):5.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。