我有以下脚本执行非线性优化（NLP），该脚本可在Matlab中运行，并在我的计算机上大约5分钟后命中MaxFunctionEvaluations：
% Generate sample consumption data (4 weeks)
x = 0:pi/8:21*pi-1e-1; %figure; plot(x, 120+5*sin(0.2*x).*exp(-2e-2*x) + 10*exp(-x))
y = 120 + 5*sin(0.2*x).*exp(-2e-2*x) + 10*exp(-x);
consumptionPerWeek = (y + [0; 11; -30; 4.5]).'; % in 168x4 format
consumptionPerHour = reshape(consumptionPerWeek, [], 1);
hoursPerWeek = 168;
hoursTotal = numel(consumptionPerHour);
daysTotal = hoursTotal/24;
weeksTotal = ceil(daysTotal/7);
%% Perform some simple calculations
q_M_mean = mean(consumptionPerHour);
dvsScalingPerWeek = mean(consumptionPerWeek)/q_M_mean;
%% Assumptions about reactor, hard-coded
V_liq = 5701.0; % m^3, main reactor; from other script
initialValue = 4.9298; % kg/m^3; from other script
substrates_FM_year = [676.5362; 451.0241];
total_DVS_year = [179.9586; 20.8867];
mean_DVS_conc = 178.1238; %kg/m^3
% Product yields (m^3 per ton DVS)
Y_M = 420;
Y_N = 389;
%% Test DVS model
DVS_hour = sum(total_DVS_year)/hoursTotal; % t/h
k_1 = 0.25; % 1/d
parameters = [k_1; Y_M; Y_N; V_liq];
%% Build reference and initial values for optimization
% Distribute feed according to demand (-24%/+26% around mean)
feedInitialMatrix = DVS_hour*ones(hoursPerWeek, 1)*dvsScalingPerWeek;
% Calculate states with reference feed (improved initials)
feedInitialVector = reshape(feedInitialMatrix, [], 1);
feedInitialVector = feedInitialVector(1:hoursTotal);
resultsRef = reactorModel1(feedInitialVector, initialValue, parameters, ...
mean_DVS_conc);
V_M_PS = 0 + cumsum(resultsRef(:,2)/24 - consumptionPerHour);
neededMStorage0 = max(V_M_PS) - min(V_M_PS);
%% Setup optimization problem (NLP): feed optimization with virtual product storage
% Objective function 1: Standard deviation of theoretical product storage volume
objFun1 = @(feedVector) objFunScalar(feedVector, initialValue, parameters, ...
mean_DVS_conc, consumptionPerHour);
% Bounds (lb <= x <= ub), i.e., decision variables can only range between 0 and 0.9*dailyDvsAmount
upperfeedLimitSlot = 0.90; % Limit DVS feed amount per *slot*
upperfeedLimitDay = 1.80; % Limit DVS feed amount per *day*
upperfeedLimitWeek = 1.37; % Limit DVS feed amount per *week*
lowerBound_nlp = zeros(1, hoursTotal);
upperBound_nlp = upperfeedLimitSlot*24*DVS_hour.*ones(1, hoursTotal);
% Equality Constraint 1: feed amount mean = constant
A_eq1_nlp = ones(1, hoursTotal);
b_eq1_nlp = DVS_hour*hoursTotal;
% Inequality Constraint 1: Limit max. daily amount
A_nlp1 = zeros(daysTotal, hoursTotal);
for dI = 1:daysTotal
A_nlp1(dI, (24*dI)-(24-1):(24*dI)) = 1;
end
b_nlp1 = upperfeedLimitDay*24*DVS_hour*ones(daysTotal, 1);
% Inequality Constraint 2: Limit max. weekly amount
A_nlp2 = zeros(weeksTotal, hoursTotal);
for wIi = 1:weeksTotal
A_nlp2(wIi, (168*wIi)-(168-1):(168*wIi)) = 1;
end
b_nlp2 = upperfeedLimitWeek*168*DVS_hour*ones(weeksTotal, 1);
% Summarize all inequality constraints
A_nlp = [A_nlp1; A_nlp2]; %sparse([A_nlp1; A_nlp2]);
b_nlp = [b_nlp1; b_nlp2]; %sparse([b_nlp1; b_nlp2]);
try
% Solver: fmincon (Matlab Optimization Toolbox) --> SQP-algorithm = best
optionen_GB = optimoptions('fmincon', 'Display', 'iter', 'FunctionTolerance', 1e-5, ...
'StepTolerance', 1e-4, 'MaxIterations', 2*hoursTotal, ...
'MaxFunctionEvaluations', 100*hoursTotal, 'HonorBounds', true, 'Algorithm', 'sqp');
catch
optionen_GB = optimset('Display', 'iter', 'TolFun', 1e-5, 'TolX', 1e-4, ...
'MaxIter', 2*hoursTotal, 'MaxFunEvals', 100*hoursTotal, 'Algorithm', 'sqp');
end
%% Solve gradient-based NLP
tic; [feedOpt, fval] = fmincon(@(feedVector) objFun1(feedVector), ...
feedInitialVector, A_nlp, b_nlp, A_eq1_nlp, b_eq1_nlp, lowerBound_nlp, upperBound_nlp, ...
[], optionen_GB); toc
%% Rerun model and calculate virtual storage volume with optimized input
resultsOpt = reactorModel1(feedOpt, initialValue, parameters, mean_DVS_conc);
q_M_Opt = resultsOpt(:,2)/24;
V_M_PS_opt = 0 + cumsum(q_M_Opt - consumptionPerHour);
neededMStorageOpt = max(V_M_PS_opt) - min(V_M_PS_opt);
sprintf('Needed product storage before optimization: %.2f m^3, \nafterwards: %.2f m^3. Reduction = %.1f %%', ...
neededMStorage0, neededMStorageOpt, (1 - neededMStorageOpt/neededMStorage0)*100)
%% Objective as separate function
function prodStorageStd = objFunScalar(dvs_feed, initialValues, parameters, mean_DVS_conc, ...
MConsumptionPerHour)
resultsAlgb = reactorModel1(dvs_feed(:, 1), initialValues, parameters, mean_DVS_conc);
q_M_prod = resultsAlgb(:,2)/24;
V_M_PS1 = 0 + cumsum(q_M_prod - MConsumptionPerHour);
prodStorageStd = std(V_M_PS1);
end
外部函数的读法如下：
function resultsArray = reactorModel1(D_feed, initialValue, parameters, D_in)
% Simulate production per hour with algebraic reactor model
% Feed is solved via a for-loop
hoursTotal = length(D_feed);
k_1 = parameters(1);
Y_M = parameters(2);
Y_N = parameters(3);
V_liq = parameters(4);
resultsArray = zeros(hoursTotal, 3);
t = 1/24;
liquid_feed = D_feed/(D_in*1e-3); % m^3/h
initialValue4Model0 = (initialValue*(V_liq - liquid_feed(1))*1e-3 ...
+ D_feed(1))*1e3/V_liq; % kg/m^3
resultsArray(1, 1) = initialValue4Model0*exp(-k_1*t);
% Simple for-loop with feed as vector per hour
for pHour = 2:hoursTotal
initialValue4Model = (resultsArray(pHour-1, 1)*(V_liq - liquid_feed(pHour))*1e-3 ...
+ D_feed(pHour))*1e3/V_liq; % kg/m^3
resultsArray(pHour, 1) = initialValue4Model*exp(-k_1*t);
end
resultsArray(:, 2) = V_liq*Y_M*k_1*resultsArray(:, 1)*1e-3; % m^3/d
resultsArray(:, 3) = V_liq*Y_N*k_1*resultsArray(:, 1)*1e-3; % m^3/d
end
当我在Octave中执行完全相同的脚本（版本5.1.0和optim 1.6.0）时，我得到：
错误：线性不平等约束：错误的尺寸
实际上是以下行（从命令提示符处执行）
sum(A_nlp*feedInitialVector <= b_nlp)
在Octave和Matlab上都给出32，因此表明尺寸正确。
这是一个错误吗？还是Octave处理线性（不等式）约束与Matlab有所不同？
（此外，如果您有如何加快此脚本的提示，它们将派上用场。）
我已对此进行了一些调试，以帮助您入门。
首先启用对错误的调试：
debug_on_error(1)
然后找到optim的安装文件夹，并查看其中的文件/private/__linear_constraint_dimensions__.m。*（我通过对您得到的确切错误进行了grep操作发现了这一点，并找到了相关的文件。在私有文件夹之外还有另一个，您可能也想看看。）] >
[如果您查看引发错误的代码行，则会注意到，例如如果rm != o.np，则触发错误，其中[rm, cm] = size(f.imc)
现在运行您的脚本，并使其在出现错误时进入调试模式。您将看到：
debug> [rm, cm] = size(f.imc) rm = 32 cm = 672 debug> o.np ans = 672 debug> rm != o.np ans = 1 % I.e. boolean test succeeds and triggers error我不知道它们是什么，大概r和c反映了行和列，但是在任何情况下，您都会发现您似乎试图将行与列进行匹配，反之亦然。
换句话说，似乎您可能在某个时候以换位方式传递了输入。
无论如何，如果这不是正在发生的事情，这应该是一个不错的起点，让您可以找出确切的错误。
我不知道为什么Matlab“有效”。也许您的代码中有一个错误，尽管它（无论好坏），matlab都能正常工作。否则，可能会出现一个错误，即无意中（或至少以与matlab不兼容的方式）对输入进行了调换。]
如果您在调试冒险后觉得它是optim软件包中的错误，请随时提交错误报告：）