suprisingly complicated optimization problem

I have the following constrained (global) optimization problem:

For a user defined sorted real values vector:

xi = [xi(1), ... , xi(N+1)]

I need to find unknown vector:

x = [x(1), ..., x(L)]

where the integer L is the unknown length of the vector x (L>=0, L = 0 is the trivial case)

and the unknown vector x = [x(1), ..., x(L)] must satisfy the following specific conditions:

===============================================================

0. The new "refined" sorted vector xo = union(xi,x), where length(xo) = L+N+1

should to fulfil the following set of conditions:

1. min(xo) = xo(1) = xi(1), max(xo) = xo(L+N+1) = xi(N+1)

2. q = max(z(j)/z(j+1),z(j+1)/z(j)) < q_max, for j = 1, 2, ... L+N

where z(j) = xo(j+1) - xo(j) , z = diff(xo), length(z) = L+N

and q_max is user defined max ratio, where q_max > 1 (typically q_max ~ 1.05 - 1.2)

3. min(z) -> maximal, minimal distance between xo vector elements should be maximized.

It is obvious that for small L the constraint conditions (2) is not possible to satisfy.

===============================================================

The motivation of this problem is the creation of the so called, "homogenized" 1-D grid, where consecutive distances between elements of vector xo are relatively "slowly" changing.

I will be very happy for any recommendation how to effectively solve this problem using MATLAB + (global) optimization toolbox.

See additional function to generate "realistic" (in my case) test input data xi:

function [x,q,x_base,x_aux] = nodegridtest(shift,show)
%
% input:
% shift ... discrete shift of aux nodes 0<=shift<=125 
% show  ... logical flag for visualization
%
% output:
% x      ... test vector
% q      ... nodes ratio vector
% x_base ... fixed part 
% x_aux  ... shifted part

% check input
if ~(0 <= shift && shift <= 125.0)
    error('shift must be in the interval [0 , 125]')
else
    shift = 2*fix(shift);
end
if nargin < 2
    show = false;
end

% realistic data generation
z_base = [7.0,7.0, ...                                    
          2.0, ...                                        
          3.4, ...                                        
          0.6, ...                                        
          7.25,7.25,7.25,7.25,7.25,7.25,7.25,7.25, ...    
          7.25,7.25,7.25,7.25,7.25,7.25,7.25,7.25, ...    
          7.25,7.25,7.25,7.25,7.25,7.25,7.25,7.25, ...    
          7.25,7.25,7.25,7.25,7.25,7.25,7.25,7.25];       
z_aux  = [3.4, ...                                       
          0.6, ...                                       
          6.0, ...                                       
          1.0, ...                                       
          8.0, ...                                       
          7.0, ...                                       
         17.5, ...                                       
          3.6, ...                                       
          6.4, ...                                       
          240.0];                                        

% upper limit of x 
xmax = sum(z_base) + sum(z_aux(1:6));   

% 
x_base        = [0, cumsum(z_base)];
x_aux         = x_base(end)+ [0, cumsum(z_aux)];
x_aux_actual  = x_aux - shift;
x_full        = union(x_base,x_aux_actual);
x             = [x_full(x_full < xmax),xmax];
x_aux         = [x_aux_actual(x_aux_actual < xmax),xmax];
z             = diff(x);
q             = max([z(1:end-1)./z(2:end);z(2:end)./z(1:end-1)],[],1);

if show
    plot(x(2:end-1),q,'bo-','LineWidth',3)
    title(['q distribution for shift = ' num2str(shift)])
    xlabel('position');ylabel('q ratio');
    hold on
    xline(x_base,'k-')
    xline(x_aux,'k--')
    hold off
    
end

then simple generate test data for any integer 0<=shift<=125 like:

xi = nodegridtest(69,1)

to get

and

xi'

ans =

         0
    7.0000
   14.0000
   16.0000
   19.4000
   20.0000
   27.2500
   34.5000
   41.7500
   49.0000
   56.2500
   63.5000
   70.7500
   78.0000
   85.2500
   92.5000
   99.7500
  107.0000
  114.0000
  114.2500
  117.4000
  118.0000
  121.5000
  124.0000
  125.0000
  128.7500
  133.0000
  136.0000
  140.0000
  143.2500
  150.5000
  157.5000
  157.7500
  161.1000
  165.0000
  167.5000
  172.2500
  179.5000
  186.7500
  194.0000
  201.2500
  208.5000
  215.7500
  223.0000
  230.2500
  237.5000
  244.7500
  252.0000
  278.0000

[Matt J]

1. min(xo) = xo(1) = xi(1), max(xo) = xo(L+N+1) = xi(N)

Do you mean max(xo) = xi(N+1)? It seems very strange that the second-last xi would set the bound.

It is obvious that for small L the constraint conditions (2) is not possible to satisfy.

Not obvious to me. If xi(N+1)/x(i)<=qmax, it seems like it should be satisfied for L=0.

[Michal]

@MattJ ad 1)see correction

Ad2) you mean xi(i+1)/x(i)<=qmax? This is trivial case, typically xi does not satisfy constraint conditions.

[Matt J]

I also don't see the difference between 2a and 2b.

[Michal]

2a and 2b are exactly identical conditions.

[Torsten]

I think the main disadvantage of all the suggested approaches is that L has to be prescribed. But I see that at least for approaches that use an optimizer turning L into a variable seems impossible because it would free the number of unknowns of the optimization.

[Matt J]

I think the main disadvantage of all the suggested approaches is that L has to be prescribed

That is only because we haven't been given a criterion on L. If the goal is to find the minimum L for which the optimization over x is feasible, then you can just use bisection, as I point out to the OP here.

[Torsten]

I'm not sure, but if you can determine the minimum L for which the corresponding grid satisfies the growth constraint, shouldn't this solution be quite near to the solution where you maximize the minimum grid spacing ?

[Michal]

Yes, proper value of L is definitely the main problem. For low L there is often no feasible solution. For high L minimum grid space maybe too low ... suprisingly complicated optimization problem.

[Michal]

OK ... after hard testing of all recommended approaches and its modifications, I can confirm that for my specific data no method is suitable to get feasible solution (if any!!!) due to the acceptable time.

Anyway, thanks to @Matt J for his effort and help.

I must continue with my search for the proper algorithm to solve this challenging optimization problem. Unknown parameter L + unknown nested vector x = [x(1), ..., x(L)] is probably the main difficulty.

Any idea how to make my optimization problem formulation more effectively solvable by available solvers?

Answer 1

[Matt J]

My thought is that you,

(1) Do not include L in your list of unknowns. Solve a series of optimization problems for fixed L=1,2,3,...

(2) Treat xo, not x, as your independent, unknown vector.

With (1) and (2), all of your listed constraints become linear, recognizing that

max(A/B,C/D)<=E

is the same as the linear constraints,

A<=B*E
C<=D*E

The only complicated constraint to express is that xi ⊂ xo. This can be done by introducing an unknown binary N+1 x (N+1+L) matrix variable P with the linear constraints,

|xi - x0|<=2*max(abs(xi)) .* (1-P)
sum(P,1)==1
sum(P,2)>=2

The objective would be rewritten as,

max q
s.t.
  q<=z  %additional linear constraints

This whole thing is an MILP, and can be solved with intlinprog. No need for a global solver. I would recommend setting it up with optimproblem.

[Michal]

xi vector must be fixed as is! This is the reason why I am looking for additional unknown vector x

[Michal]

Could you more elaborate your suggestion? Condition (2) does not look as linear.

max(A/B,C/D)<=E

condition (2)should looks like:

max(A/B,B/A)<=E

What do you mean by q in your objective?

[Matt J]

xi vector must be fixed as is! This is the reason why I am looking for additional unknown vector x

xi is still fixed in my formulation. It is just that none of the unknowns xo(k) are fixed.

condition (2)should looks like: max(A/B,B/A)<=E

That's the same as what I wrote except you've made the substitutions C=B, D=A. It leads to the linear constraints,

A<=B*E
B<=A*E

What do you mean by q in your objective?

It's an additional unknown variable. So your total set of optimization variables will be (q,P, xo)

[Michal]

@Matt J I still don't fully understand your suggestions. They are too general. Could you show me how to do this with a suitable specific solver?

[Michal]

Do not include L in your list of unknowns. Solve a series of optimization problems for fixed L=1,2,3,...

OK, this could be a simpler way. Problem is, that for some xi vector could be L>>1

[Matt J]

OK, this could be a simpler way. Problem is, that for some xi vector could be L>>1

Isn't the optimal min(z) a decreasing function of L? Are we not just trying to find the smallest L for which the optimization is feasible? If so, you could use bi-section on L rather than a linear search L=1,2,3,...

[Michal]

Ones again, could you briefly show me how to implement proper solver on this task (L=1,2,3,...). I am not so familiar with Optimization toolbox solvers.

[William Rose]

Very impressive answer by @Matt J

[Michal]

I am still confused, because solver recommended by @Matt J intlinprog does not offer nonlinear constraints which looks like:

|xi - xo|<=2*max(abs(xi)) .* (1-P)

Of course, could be possible that I again overlooked something important :) How to rewrite this constraints as A·x ≤ b?

Is there anybody willing to show me how to implement this optimization problem with the proper solver?

[Matt J]

Here's an example:

%Problem data
xi=sort(linspace(0,1,5)');
L=10;
qmax=2;

N=numel(xi)-1;
U=max(xi);

eNL=ones(N+L+1,1);
eN=ones(N+1,1);


%Unknowns
q =  optimvar('q');
P =  optimvar('P',N+1, L+N+1, Type='integer', Lower=0,Upper=1);
xo = optimvar('xo',L+N+1, Lower=0,Upper=U);

z=diff(xo);

%Constraints
Constraints.con_mono=z>=0; %monotonicity of xo
Constraints.conq= q<=z;   %defines object value q 

Constraints.con2a_1=z(1:end-1)<=qmax*z(2:end);  %qmax constraints
Constraints.con2a_2=z(2:end)<=qmax*z(1:end-1);

Constraints.conP_1=xi*eNL'- eN*xo' <=2*U*(1-P); %constraints on P
Constraints.conP_2=eN*xo' - xi*eNL'<=2*U*(1-P);
Constraints.conP_3=sum(P,2)==1;
Constraints.conP_4=sum(P,1)<=1;

%Solve
prob=optimproblem('Objective',q,'ObjectiveSense','max','Constraints', Constraints);

sol = solve(prob)

Solving problem using intlinprog.
Running HiGHS 1.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 3e+00]
  Cost   [1e+00, 1e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 3e+00]
Presolving model
194 rows, 91 cols, 538 nonzeros  0s
188 rows, 47 cols, 374 nonzeros  0s
132 rows, 47 cols, 424 nonzeros  0s

Solving MIP model with:
   132 rows
   47 cols (33 binary, 0 integer, 0 implied int., 14 continuous)
   424 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0.6666666667    -inf                 inf        0      0      0         0     0.0s
         0       0         0   0.00%   0.0714285714    -inf                 inf        0      0      5        25     0.0s
 R       0       0         0   0.00%   0.0714285714    0.05              42.86%      115      6     20        39     0.0s
 L       0       0         0   0.00%   0.0714285714    0.0625            14.29%      117      7     24        46     0.0s

48.5% inactive integer columns, restarting
Model after restart has 75 rows, 24 cols (8 bin., 0 int., 0 impl., 14 cont.), and 188 nonzeros

         0       0         0   0.00%   0.0714285714    0.0625            14.29%        4      0      0        64     0.0s
         0       0         0   0.00%   0.0714285714    0.0625            14.29%        4      0      0        69     0.0s

Solving report
  Status            Optimal
  Primal bound      0.0625
  Dual bound        0.0625
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0.0625 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.03 (total)
                    0.01 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     69 (total)
                    0 (strong br.)
                    28 (separation)
                    11 (heuristics)

Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are
integer within tolerance, options.ConstraintTolerance = 1e-06.

sol = struct with fields:
     P: [5×15 double]
     q: 0.0625
    xo: [15×1 double]

xo=sol.xo'

xo = 1×15

0 0.0625 0.1250 0.2500 0.3125 0.3750 0.5000 0.5625 0.6250 0.6875 0.7500 0.8125 0.8750 0.9375 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xi'

ans = 1×5

0 0.2500 0.5000 0.7500 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Michal]

Thanks a lot @Matt J! Now is situation much more clear...

There are two problems:

1. condition (1) is not satisfied reliably when I create the data differently

xi=[sort(rand(1,5))]';
L=10;
qmax=1.2;

with output:

>> xi'

ans =

0.4501    0.4642    0.4816    0.5246    0.5315
>> sol.xo'

ans =

0.4033    0.4183    0.4309    0.4413    0.4501    0.4573    0.4642    0.4725    0.4816    0.4892    0.4970    0.5063    0.5163    0.5246    0.5315

you can see first 4 wrong elements of sol.xo

But the 2nd problem is more serious:

2. for typical value is N ~ 50-100 and L ~ (2-5)*N is HiGHS algorithm in intlinprog function terribly slow. Probably due to the size of binary marix P (N+1) x (N+1+L) and correponding constraints and/or complexity of the used MILP algorithm.

Any idea how to speed up solver?

[Matt J]

1. condition (1) is not satisfied reliably when I create the data differently

It's because xo.Lower was hardcoded to 0. I've fixed it in the revision below.

Any idea how to speed up solver?

Not sure. Maybe by providing initial guesses for the unknowns? The example below includes one way of doing this, but you would have to try it out on realistic data sizes.

%Problem data
xi=[sort(rand(1,5))]';
L=10;
qmax=1.2;

N=numel(xi)-1;
xiMax=max(xi);
xiMin=min(xi);
U=max(abs(xiMin), abs(xiMax));

eNL=ones(N+L+1,1);
eN=ones(N+1,1);


%Unknowns
q =  optimvar('q');
P =  optimvar('P',N+1, L+N+1, Type='integer', Lower=0,Upper=1);
xo = optimvar('xo',L+N+1, Lower=xiMin,Upper=xiMax);

z=diff(xo);

%Constraints
Constraints.con_mono=z>=0; %monotonicity of xo
Constraints.conq= q<=z;    %defines objective value q 

Constraints.con2a_1=z(1:end-1)<=qmax*z(2:end);  %qmax constraints
Constraints.con2a_2=z(2:end)<=qmax*z(1:end-1);

Constraints.conP_1=xi*eNL'- eN*xo' <=2*U*(1-P); %constraints on P
Constraints.conP_2=eN*xo' - xi*eNL'<=2*U*(1-P);
Constraints.conP_3=sum(P,2)==1;
Constraints.conP_4=sum(P,1)<=1;

%Solve
prob=optimproblem('Objective',q,'ObjectiveSense','max','Constraints', Constraints);

%Possible initial guess
y=linspace(xiMin,xiMax,L+2); y([1,end])=[];
sol0.xo=union(xi',y)';
sol0.P=double(xi==sol0.xo');
sol0.q=min(diff(sol0.xo));

sol = solve(prob,sol0)

Solving problem using intlinprog.
Running HiGHS 1.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [1e+00, 1e+00]
  Bound  [7e-02, 1e+00]
  RHS    [5e-01, 2e+00]
Assessing feasibility of MIP using primal feasibility and integrality tolerance of       1e-06
WARNING: Row     33 has         infeasibility of    0.002658 from [lower, value, upper] = [           -inf;    0.0026583052;               0]
WARNING: Row     34 has         infeasibility of     0.02126 from [lower, value, upper] = [           -inf;     0.021261431;               0]
Solution has               num          max          sum
Col     infeasibilities      0            0            0
Integer infeasibilities      0            0            0
Row     infeasibilities      8      0.03556       0.1766
Row     residuals            0            0            0
Attempting to find feasible solution for (partial) user-supplied values of discrete variables
Coefficient ranges:
  Matrix [1e+00, 2e+00]
  Cost   [1e+00, 1e+00]
  Bound  [7e-02, 1e+00]
  RHS    [5e-01, 2e+00]
Presolving model
58 rows, 14 cols, 144 nonzeros  0s
44 rows, 11 cols, 100 nonzeros  0s
Presolve : Reductions: rows 44(-180); columns 11(-80); elements 100(-498)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.0000000000e+03 Ph1: 14(14000); Du: 1(1) 0s
         18    -2.9095257068e-02 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 18
Objective value     :  2.9095257068e-02
HiGHS run time      :          0.00
Presolving model
194 rows, 91 cols, 538 nonzeros  0s
186 rows, 41 cols, 348 nonzeros  0s
114 rows, 33 cols, 270 nonzeros  0s
94 rows, 32 cols, 246 nonzeros  0s

MIP start solution is feasible, objective value is 0.0290952570677

Solving MIP model with:
   94 rows
   32 cols (18 binary, 0 integer, 0 implied int., 14 continuous)
   246 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0.1964387286    0.0290952571     575.16%        0      0      0         0     0.0s
         0       0         0   0.00%   0.0324231086    0.0290952571      11.44%        0      0      4        22     0.0s

Solving report
  Status            Optimal
  Primal bound      0.0290952570677
  Dual bound        0.0290952570677
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0.0290952570677 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.01 (total)
                    0.01 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     51 (total)
                    0 (strong br.)
                    9 (separation)
                    0 (heuristics)

Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are
integer within tolerance, options.ConstraintTolerance = 1e-06.

sol = struct with fields:
     P: [5×15 double]
     q: 0.0291
    xo: [15×1 double]

xo=sol.xo'

xo = 1×15

0.0677 0.1062 0.1411 0.1702 0.1993 0.2323 0.2720 0.3069 0.3360 0.3651 0.3994 0.4343 0.4634 0.4925 0.5216

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xi'

ans = 1×5

0.0677 0.1411 0.3069 0.4343 0.5216

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Michal]

thanks for fix ...

regarding speed up by providing the initial guess as:

%Possible initial guess
y=linspace(xiMin,xiMax,L+2); y([1,end])=[];
sol0.xo=union(xi',y)';
sol0.P=double(xi==sol0.xo');
sol0.q=min(diff(sol0.xo));

is not viable way ... no significant speed up :(

Answer 2

[Matt J]

Here's another formulation, closer to your original one. It uses fmincon and is really fast (when the problem is feasible), but one of the constraints is nonconvex and only piecewise differentiable. You will probably have to incorporate MultiStart to help ensure global optimality.

%Problem data
xi=[sort(rand(1,5))]';
L=10;
qmax=1.2;

N=numel(xi)-1;
xiMax=max(xi);
xiMin=min(xi);
U=max(abs(xiMin), abs(xiMax));

eNL=ones(N+L+1,1);
eN=ones(N+1,1);


%Unknowns
x = optimvar('x', L, Lower=xiMin,Upper=xiMax);
z = optimvar('z', L+N, Lower=0,Upper=xiMax-xiMin);
fobj = optimvar('fobj', Lower=0);


%Constraints
fcn=@(c)diff(sort([c(:);xi(:)]));
Constraints.con_zdef= z(:) == fcn2optimexpr(fcn,x);

Constraints.con2a_1=z(1:end-1)<=qmax*z(2:end);  %qmax constraints
Constraints.con2a_2=z(2:end)<=qmax*z(1:end-1);

Constraints.obj=fobj<=z;


%Solve
prob=optimproblem('Objective',fobj,'ObjectiveSense','max','Constraints', Constraints);

sol0.x=linspace(xiMin,xiMax,L)';
sol0.z=fcn(sol0.x);
sol0.fobj=min(sol0.z);

opts=optimoptions('fmincon','StepTol',1e-10,'OptimalityTol',1e-10,'FunctionTol',1e-10, ...
                        'MaxFunEvals',1e7,'MaxIterations',1e5);

sol = solve(prob, sol0,Options=opts)

Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol = struct with fields:
    fobj: 0.0238
       x: [10×1 double]
       z: [14×1 double]

xo=union(xi,sol.x)'

xo = 1×15

0.3218 0.3657 0.4084 0.4563 0.5120 0.5757 0.6460 0.7147 0.7759 0.8292 0.8750 0.9141 0.9473 0.9754 0.9992

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xi'

ans = 1×5

0.3218 0.3657 0.4563 0.9754 0.9992

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Michal]

I do not know what is wrong but this code is terribly slow. Probably due to the too tight parameters:

opts=optimoptions('fmincon','StepTol',1e-10,'OptimalityTol',1e-10,'FunctionTol',1e-10,'MaxFunEvals',1e7,'MaxIterations',1e5);

[Matt J]

No, the output above took a blink of an eye for me. The problem is that sometimes the randomly generated problem data doesn't give an optimization problem with a feasible solution. It can take fmincon a very long time to figure out that the optimization problem is infeasible.

[Matt J]

It can take fmincon a very long time to figure out that the optimization problem is infeasible.

You can abort the optimization if it runs to long using an OuputFcn,

https://www.mathworks.com/help/optim/ug/output-function.html

Answer 3

[Matt J]

This is another suggested strategy: to find the minimum possible qmax, given L. This seems at the very least like an important first step, since you need to see if the qmax values you're pursuing are even achievable.

I demonstrate using the xi data generated in your posted example,

xi = nodegridtest(69,1)

It uses minL1intlin, as suggested earlier, to derive an initial guess.

clear,clc,close all

%Problem data
load xiData

tic;
[~,qmax]=runOptimization(xi, numel(xi))

Warning: Cannot use sparse matrices with sqp algorithm: converting to full.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

qmax = 4.8151

toc;

Elapsed time is 3.221954 seconds.

tic;
[~,qmax]=runOptimization(xi, 5*numel(xi))

Warning: Cannot use sparse matrices with sqp algorithm: converting to full.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

qmax = 1.7836

toc;

Elapsed time is 0.857964 seconds.

function [xo,qmax]=runOptimization(xi,L)
% Inputs:
%     xi - Description of input xi
%     L - Description of input L
%
% Outputs:
%     xo - Description of output xo
%     qmax - Description of output qmax
arguments
 xi   (:,1) double
 L    (1,1) double {mustBeInteger}
end

    N=numel(xi)-1;
    xiMax=max(xi);
    xiMin=min(xi);

    fcn=@(c)log(diff(sort([c(:);xi(:)])));

    %Unknowns
    x = optimvar('x', L, Lower=xiMin,Upper=xiMax);
    logz = optimvar('logz', L+N, Lower=-inf,Upper=log(xiMax-xiMin));
    y = optimvar('y', Lower=0); %log(qmax)


    %Constraints

    Constraints.qmax_1=diff(logz)<=y;  %qmax constraints
    Constraints.qmax_2=-diff(logz)<=y;

    Constraints.nlcon=x(1)^2<=1; %dummy

    %Solve
    prob=optimproblem('Objective',y,'ObjectiveSense','min','Constraints', Constraints);



    sol0.x=initialX(L,xi);
    sol0.logz=fcn(sol0.x);
    sol0.y=norm(sol0.logz,inf);




    opts=optimoptions('fmincon','StepTol',1e-6,'OptimalityTol',1e-6,'FunctionTol',1e-6, ...
                            'MaxFunEvals',1e7,'MaxIterations',1e3,'SpecifyObjectiveGradient',true, ...
                            'SpecifyConstraintGradient',true,'Display','final','Algorithm','sqp');

    Prob=prob2struct(prob,sol0,'Solver','fmincon');
    Prob.nonlcon=@(vars)nlcon(vars,xi,N,L); %
    Prob.options=opts;


    vars = fmincon(Prob);


       m=N+L+1;
       vars=vars(:);
       logz=vars(1:m-1); %#ok<NASGU>
       x=vars(m+1:m+L);
       y=vars(m);

       qmax=exp(y);

    xo=union(xi,x)';

end


function [c,ceq,cgrad, ceqgrad] = nlcon(vars,xi,N,L)

   c=[]; cgrad=[];

   m=N+L+1;

   vars=vars(:);
   logz=vars(1:m-1);
   x=vars(m+1:m+L);

   [xo,is]=sort([x;xi]);
   z=diff(xo);

   ceq = logz - log(z);

   if nargout>3

       P=sparse(1:m,is,1, m,m);
       P=P(:,1:L);

       Jxo=diff(P,1,1)./z;

       ceqgrad=[speye(m-1,m), -Jxo].';

   end

end

function x=initialX(L,xi)

 N=numel(xi)-1;
 xiMax=max(xi);
 xiMin=min(xi);

 dL0 = diff(xi)./(xiMax-xiMin)*L;
 opts=optimoptions('intlinprog','Display','none');
 dL=minL1intlin(speye(N),round(dL0(:)),1:N,[],[],ones(1,N),L,zeros(N,1),...
                [],[],opts );

 K=numel(xi)-1;
 xo=cell(1,K);
 for k=1:numel(xi)-1
  tmp=linspace(xi(k), xi(k+1),dL(k)+2);
  xo{k}=tmp(1:end-1);
 end
 xo=[cell2mat(xo),xi(end)];

 x=setdiff(xo,xi);


end

[Michal]

@Matt J Thanks...!!! I need to test your code more thoroughly.

Answer 4

[Matt J]

The motivation of this problem is the creation of the so called, "homogenized" 1-D grid, where consecutive distances between elements of vector xo are relatively "slowly" changing.

If so, you might consider a simplified solution using minL1intlin from this File Exchange submission,

https://www.mathworks.com/matlabcentral/fileexchange/52795-constrained-minimum-l1-norm-solutions-of-linear-equations?s_tid=srchtitle

Basically, this chooses a number of points dL(k) to place between each xi, with sum(dL)==L and in proportion to the spacing xi(k+1)-xi(k). I think this is almost what you want, though there is no qmax parameter to let you bound the z ratios. If nothing else, it may give you a good way to form an initial guess for one of the other solution methods.

%Problem data
xi=[sort(rand(1,7))]';
L=12;
% qmax=1.2;

N=numel(xi)-1;
xiMax=max(xi);
xiMin=min(xi);

dL0 = diff(xi)./(xiMax-xiMin)*L;

dL=minL1intlin(speye(N),round(dL0(:)),1:N,[],[],ones(1,N),L,zeros(N,1) );

Running HiGHS 1.7.1: Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 1e+00]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+01]
Presolving model
9 rows, 10 cols, 22 nonzeros  0s
9 rows, 9 cols, 21 nonzeros  0s
Objective function is integral with scale 1

Solving MIP model with:
   9 rows
   9 cols (0 binary, 5 integer, 4 implied int., 0 continuous)
   21 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s
 T       0       0         0   0.00%   0               1                100.00%        0      0      0         7     0.0s

Solving report
  Status            Optimal
  Primal bound      1
  Dual bound        1
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    1 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            0.01 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             1
  LP iterations     7 (total)
                    0 (strong br.)
                    0 (separation)
                    0 (heuristics)

Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are
integer within tolerance, options.ConstraintTolerance = 1e-06.

xo=getxo(dL,xi);

check = numel(xo) == L+N+1

check = logical
   1

h=plot(xi,0*xi,'x',xo,0*xo,'o'); axis padded
set(h,'MarkerSize',8);
legend xi xo

●

function xo=getxo(dL,xi)

 K=numel(xi)-1;
 xo=cell(1,K);
 for k=1:numel(xi)-1

  tmp=linspace(xi(k), xi(k+1),dL(k)+2);
  xo{k}=tmp(1:end-1);

 end

 xo=[cell2mat(xo),xi(end)];
  

end

[Michal]

I think that some control parameters (qmax, zmin, ...) for acceptable grid inhomogenity are always very important.

[Matt J]

But in your formulation, the user is not free to set these parameters because it is obscure which combinations of parameters give a feasible problem. The simplified formulation above is always feasible.

[Michal]

Yes, you right. But this solver produce nearly always unacceptable high ratios q > qmax. And qmax is really important for me!

[Matt J]

And qmax is really important for me!

As shown below, (with your posted xi data, and L=3*N) even though the max q ratio is high with this strategy, the majority of the ratios are very low, and the distribution of spacings z(i) has a low variance. Are you sure that's not enough?

%Problem data
load xiData


N=numel(xi)-1;

L=3*N;

xiMax=max(xi);
xiMin=min(xi);

dL0 = diff(xi)./(xiMax-xiMin)*L;

opts=optimoptions('intlinprog','Display','none');
dL=minL1intlin(speye(N),round(dL0(:)),1:N,[],[],ones(1,N),L,zeros(N,1),[],[],opts );

xo=getxo(dL,xi);

z=diff(xo);
q=z(2:end)./z(1:end-1);
q=max(q,1./q);

tiledlayout('h');
nexttile; histogram(q,100); title q-Histogram; axis square
nexttile; histogram(z,100); title z-Histogram; axis square

function xo=getxo(dL,xi)

 K=numel(xi)-1;
 xo=cell(1,K);
 for k=1:numel(xi)-1

  tmp=linspace(xi(k), xi(k+1),dL(k)+2);
  xo{k}=tmp(1:end-1);

 end

 xo=[cell2mat(xo),xi(end)];
  

end