3

I have the following script performing a nonlinear optimization (NLP), which works in Matlab and hits MaxFunctionEvaluations after about 5 minutes on my machine:

% 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

The external function reads like this:

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

When I execute the very same script in Octave (ver 5.1.0 with optim 1.6.0), I get:

error: linear inequality constraints: wrong dimensions

When in fact, the following line (executed from the command prompt)

sum(A_nlp*feedInitialVector <= b_nlp)

gives 32 on both Octave and Matlab, thus showing that dimensions are correct.

Is this a bug? Or is Octave treating linear (in)equality constraints somehow different than Matlab?

(Also, if you have tips on how to speed up this script, they would come in handy.)

winkmal
  • 622
  • 1
  • 6
  • 16
  • Actually, I think this is a bug. Did you have more rows than columns in your constraints? – Jan Oct 15 '20 at 11:54
  • If I understand you correctly, the answer is no. `A_nlp` is of shape `(32x672)`, and `A_eq1_nlp` `(1x672)`. So there are much more columns than rows. – winkmal Oct 15 '20 at 13:05
  • OK -- my matrix is like `(2000x10)`. So this is not the problem. Still I find this very strange since in Matlab it works just fine. – Jan Oct 15 '20 at 14:30

1 Answers1

3

I've debugged this a bit for you to get you started.

First enable debugging on error:

debug_on_error(1)

Then find the installation folder of optim, and have a look at file /private/__linear_constraint_dimensions__.m within.
*(I found this by doing a grep operation for the exact error you were getting, and found the relevant file. There is another one outside the private folder, you may want to look at that too.)

If you look at the lines trigerring the errors, you will notice, e.g. that an error is triggered if rm != o.np, where [rm, cm] = size(f.imc)

Now run your script and let it enter debug mode on error. You will see that:

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

I have no idea what these are, presumably r and c reflect rows and columns, but in any case, you will see that it appears you are trying to match rows with columns and vice versa.

In other words, it looks like you may have passed your inputs in a transposed fashion at some point.

In any case, if this isn't exactly what's happening, this should be a decent starting point for you to figure the exact bug out.

I don't know why matlab "works". Maybe there's a bug in your code and matlab works despite it (for better or worse).
Or there might be a bug in optim transposing inputs by accident (or, at least, in a manner that is incompatible to matlab).

If you feel after your debugging adventures that it's a bug in the optim package, feel free to file a bug report :)

Tasos Papastylianou
  • 21,371
  • 2
  • 28
  • 57
  • 2
    I actually got it to work by *transposing* `A_*`. I do not find this intuitive at all, since `A_nlp*feedInitialVector` works, but `A_nlp.' * feedInitialVector` throws an error. Also, I had to move the objective function from the end of the script to the middle, or a separate function. Finally, `fmincon` does not start since **Initial parameters violate constraints**. The difference between left- and right-hand side `A_eq1_nlp*feedInitialVector - b_eq1_nlp` is `1.4211e-13`... If I manually add this to `b_eq1_nlp`, it finally starts to run, though it takes like 18x to solve than in Matlab. – winkmal Jan 30 '20 at 13:49
  • @rotton glad you worked it out. For me the main concern is the last one you mentioned, i.e. speed. Unfortunately, that is one area where Octave is still behind matlab, but I suppose this lag is the temporary price to pay for having an otherwise high quality free and open-source alternative to a proprietary product. :) Regarding the other claims though, I would say that these are arguments *for* octave rather than against. If it is stricter, and less forgiving in terms of expecting arguments in proper form and strictly adhering to constraints, then this is good, and will save you from errors! – Tasos Papastylianou Jan 30 '20 at 15:37
  • @rotton the point about having to move the function further up, this is one of the big (in my view reasonable) differences between octave and matlab. In octave, a _local_ function needs to have been defined before use. This allows you to redefine a function safely. In other words, octave is more like python and julia in this regard. Matlab has recently changed the behaviour of local functions, such that they can be defined anywhere in the file, or explicitly at the end in the case of scripts. In any case, the middle ground is better style anyway: use subfunctions, not locally defined ones. – Tasos Papastylianou Jan 30 '20 at 15:41