为什么这些线性不等式约束在 Matlab 中有效但在 Octave 中无效？
Why do these linear inequality constraints work in Matlab but not Octave?
我有以下脚本执行非线性优化 (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）中执行完全相同的脚本时，我得到：
error: linear inequality constraints: wrong dimensions
实际上，以下行（从命令提示符执行）
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 "works"。也许您的代码中存在错误，尽管如此，matlab 仍能正常工作（无论好坏）。
或者在 optim 转置输入中可能存在错误（或者，至少以与 matlab 不兼容的方式）。
如果您在调试之后觉得这是 optim 包中的错误，请随时提交错误报告:)
我有以下脚本执行非线性优化 (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）中执行完全相同的脚本时，我得到：
error: linear inequality constraints: wrong dimensions
实际上，以下行（从命令提示符执行）
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 "works"。也许您的代码中存在错误，尽管如此，matlab 仍能正常工作（无论好坏）。
或者在 optim 转置输入中可能存在错误（或者，至少以与 matlab 不兼容的方式）。
如果您在调试之后觉得这是 optim 包中的错误，请随时提交错误报告:)