У меня есть следующий скрипт, выполняющий нелинейную оптимизацию (NLP), который работает в Matlab и нажимает MaxFunctionEvaluations примерно через 5 минут на моей машине:
% 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)
дает 32 как для Octave, так и для Matlab, показывая, таким образом, правильные размеры .
Это ошибка? Или Octave рассматривает линейные (не) ограничения равенства как-то иначе, чем Matlab?
(Кроме того, если у вас есть советы о том, как ускорить этот сценарий, они пригодятся.)