Numeric - symbolic solution of GPS phase ambiguity problem with Mathematica

Béla Paláncz

January 9, 2018

© Béla Paláncz 2018

This contribution consists of four parts. All parts contains the Mathematica codes as well as their explanations. The codes can be generally employed for other examples, too. For this demonstration Mathematica Version 11.2 running under MS Windows 7 on computer HP Z420 has been employed.

Contents

Abstract

1. Linearization of the problem

1.1 Brute force solution

1.2 Computing the bounds

1.3 McCormic envelopes

1.3.1 Convex linear relaxation

1.3.2 Application of McCormic envelopes

2. Solving linear problem

2.1 Nonlinear methods

2.1.1 Integer nonlinear optimization

2.1.2 Global numerical method

2.1.3 Local numerical method

2.2 Linear programming

2.3 Successive approximation

3. Mixed integer problem

3.1 Converting mixed problem into pure one

4. GPS phase ambiguity problem

4.1 Problem description

4.2 Input data

4.3 Data preprocessing

4.4 Solution methods

4.4.1 Direct nonlinear programming

4.4.2 Continuous solution

4.4.3 A naive solution

4.4.4 Linearized problem

4.4.5 Iterative solution

4.4.6 Solution without iteration

4.4.7 Employing linear programming

5. Conclusions

Literature

The GPS phase ambiguity is a global, quadratic mixed integer programming problem, which should be computed online. To solve the problem, probably the most popular procedure is the so called LAMBDA method (see Teunissen et al, 1997). In this article we suggest an alternative algorithm, which can also satisfy this requirement. The algorithm utilizes the Mathematica’s numeric - symbolic computation ability, first transforms the mixed problem into a pure integer one. Then the bounds of region of the global solution can be computed on the basis of maximum - minimum eigenvalues of the matrix of the bilinear form with the continuous solution as center of this region. Employing McCormic Envelopes the quadratic problem is linearized. As illustration of the method a simple example as well as a real world five - satellite problem are given.

Contents

1.1 Brute force solution

1.2 Computing the bounds

1.3 McCormic envelopes

1.3.1 Convex linear relaxation

1.3.2 Application of McCormic envelopes

1.1 Brute force solution

In order to make the algorithm more transparent, we demonstrate the steps of the solution of a simple global quadratic integer programming problem.

Let us consider the following example adapted from Li and Sun (2006). We are looking for the global integer minimum of this objective,

q=27 ((Subscript[x, 1])^2)-18 (Subscript[x, 1]) (Subscript[x, 2])+4 ((Subscript[x, 2])^2)-3 (Subscript[x, 2]);

The contour plot of the objective is,

p0=ContourPlot[q,{(Subscript[x, 1]),-4,6},{(Subscript[x, 2]),-4,6},Frame->True,FrameLabel ->Automatic,ContourShading->None,MeshStyle->{Red,Blue},MeshFunctions->{ #1&,#2&},Mesh->{9,9},ContourLabels->True,AspectRatio->1,PlotRange->{{-4,6},{-4,6}}]

Fig .1 The contour plot of the objective

As a first approach we try to solve the  problem via nonlinear integer programming. Excluding the trivial {0,0} solution,

Minimize[{q,(Subscript[x, 1])≠0,Element[{(Subscript[x, 1]),(Subscript[x, 2])},Integers]},{(Subscript[x, 1]),(Subscript[x, 2])}]//AbsoluteTiming

It looks an easy way, however in general, the global, unconstrained mixed quadratic integer programming problem is an NP - hard problem and it is not possible to solve directly in practical cases, when problem size is large. Therefore our aim is to reduce the problem to a linear integer programming task. Since linear optimization needs constraints for the variables, the first step is to compute the bounds of the global solution.

1.2 Computing the bounds

These bounds can be computed using  the continuous global solution.The optimal solution is,

NMinimize[q,{(Subscript[x, 1]),(Subscript[x, 2])}]//AbsoluteTiming

Let us consider the bilinear form of the objective with the continuous solution as center in the form,

(({{(Overscript[(Subscript[x, 1]), ∼])-(1/2), (Overscript[(Subscript[x, 2]), ∼])-(3/2)}}))𝒬(({{(Overscript[(Subscript[x, 1]), ∼])-(1/2)}, {(Overscript[(Subscript[x, 2]), ∼])-(3/2)}}))

where (Overscript[(Subscript[x, 1]), ∼]) and (Overscript[(Subscript[x, 2]), ∼]) are the integer solutions of the optimization problem. The matrix 𝒬 can be computed as

𝒬=(1/2)(({{(((∂^2)q)/(∂((Subscript[x, 1])^2))), (((∂^2)q)/(∂(Subscript[x, 1])∂(Subscript[x, 2])))}, {(((∂^2)q)/(∂(Subscript[x, 1])∂(Subscript[x, 2]))), (((∂^2)q)/(∂((Subscript[x, 2])^2)))}}))

Then

Q=1/2 D[q,{{(Subscript[x, 1]),(Subscript[x, 2])},2}];Q//MatrixForm

Now computing  the eigenvalues,

λ=Eigenvalues[Q]//N

The ratio of the maximum  and minimum values of λ’s

κ=((Max[λ])/(Min[λ]))

Then the radius of a n - dimensional hyper sphere with the continuous solution as center, R = 1/2 (Sqrt[n κ,]) now n = 2,

R=1/2 (Sqrt[2 κ])

Let us illustrate the bound of the global solution. The green point represents the continuous solution.

p1=Show[{p0,ListPlot[{{1/2,3/2}},PlotStyle->{Green,PointSize[0.02]},AspectRatio->1],ContourPlot[(R^2)==((Subscript[x, 1])-1/2)^2+((Subscript[x, 2])-3/2)^2,{(Subscript[x, 1]),-4,6},{(Subscript[x, 2]),-4,6},AspectRatio->1]}]

Fig. 2 The bounded region of the global integer optimum

We are looking for a box - type constrain, then the coordinates of the lower and upper bounds can be computed as

(((u-1/2))^2)+(((v-3/2))^2)-(R^2)/.v->3/2

Solve[%==0,u]

and

(((u-1/2))^2)+(((v-3/2))^2)-(R^2)/.u->1/2

Solve[%==0,v]

Then the box constraint,

-3≤(Subscript[x, 1])≤4,-2≤(Subscript[x, 2])≤5

Show[{p1,Graphics[{Yellow,Rectangle[{-3,-2},{4,5}]},Frame->True],p1}]

Fig .3 The box - bounded region of the global integer optimum

Now, we try to solve the constrained quadratic integer problem directly via nonlinear programming. The constrains avoiding {0,0} trivial solution,

constraints=Apply[And,{ -3≤(Subscript[x, 1])≤4,-2≤(Subscript[x, 2])≤5,(Subscript[x, 1])≠0,(Subscript[x, 2])≠0,Element[{(Subscript[x, 1]),(Subscript[x, 2])},Integers]}]

Minimize[{q,constraints},{(Subscript[x, 1]),(Subscript[x, 2])}]//AbsoluteTiming

Having the bounds, the computation time is considerably shorter, however our aim is to linearize the quadratic objective. We shall employ the McCormick Envelopes.

1.3 McCormic envelopes

The McCormick envelopes the convex relaxation of a quadratic problem via introducing new variables for the quadratic terms and employing the  additional constraints. In general, we introduce new variables,

(Subscript[w , ij ])=(Subscript[x, i])(Subscript[x, j])

with following constraints,

(Subscript[w, ij])≥(Subscript[x, iL]) (Subscript[x, j])+(Subscript[x, i])(Subscript[x, jL])-(Subscript[x, iL]) (Subscript[x, jL]),
(Subscript[w, ij])≥(Subscript[x, iU]) (Subscript[x, j])+(Subscript[x, i])(Subscript[x, jU])-(Subscript[x, iU]) (Subscript[x, jU]),
(Subscript[w, ij])≤(Subscript[x, iU]) (Subscript[x, j])+(Subscript[x, jL])(Subscript[x, i])-(Subscript[x, iU]) (Subscript[x, jL]),
(Subscript[w, ij])≤(Subscript[x, i])(Subscript[x, jU])+(Subscript[x, iL]) (Subscript[x, j])-(Subscript[x, jU]) (Subscript[x, iL]),

where (Subscript[x, L])≤ x ≤(Subscript[x, U]).

1.3.1 Convex linear relaxation

Let us illustrate this technique with a very simple 1D example. The quadratic objective is  y(x) = (x^2), where the constraints are  -1 ≤ x ≤ 2. Then introducing  w = (x^2)

(Subscript[x, L])=-1;(Subscript[x, U])=2;

and the additional constraints,

{c1,c2,c3,c4}={w≥(Subscript[x, L]) x+x (Subscript[x, L])-(Subscript[x, L]) (Subscript[x, L]),w≥(Subscript[x, U]) x+x (Subscript[x, U])-(Subscript[x, U]) (Subscript[x, U]),w≤(Subscript[x, U]) x+x (Subscript[x, L])-(Subscript[x, U]) (Subscript[x, L]),(Subscript[x, L])≤x≤(Subscript[x, U])}

we have a linear system. Let solve it,

solx=Minimize[{w,c1,c2,c3,c4},{w,x}]

Displaying the linear constraints (lines) it is clear how this convex linear envelopes “embrace” the objective and provide and approximate solution of the quadratic problem.

Show[{Plot[{(x^2),c1[[2]],c2[[2]],c3[[2]]},{x,-1,2},AspectRatio->1.5,PlotRange->{-3,4},AxesLabel->{x,w}],ListPlot[{{solx[[2,2,2]],solx[[1]]}},PlotStyle->{PointSize[0.03],Red}]}]

Fig .4 The McCormick envelopes (triangular region) and the minimum of the linear problem (red point)

This is the first linear approximation of the solution. A better approximation can be achieved when  we reduce the constraints  using the solution ( x = (1/2)), as new upper limit,

(Subscript[x, L])=-1;(Subscript[x, U])=(1/2);

and the new linear constrains (envelopes) are,

{c1,c2,c3,c4}={w≥(Subscript[x, L]) x+x (Subscript[x, L])-(Subscript[x, L]) (Subscript[x, L]),w≥(Subscript[x, U]) x+x (Subscript[x, U])-(Subscript[x, U]) (Subscript[x, U]),w≤(Subscript[x, U]) x+x (Subscript[x, L])-(Subscript[x, U]) (Subscript[x, L]),(Subscript[x, L])≤x≤(Subscript[x, U])}

Then the result of the next approach,

solx=Minimize[{w,c1,c2,c3,c4},{w,x}]

So we have got closer to x = 0.

Show[{Plot[{(x^2),c1[[2]],c2[[2]],c3[[2]]},{x,-1,(1/2)},AspectRatio->1.5,PlotRange->{-0.75,1},AxesLabel->{x,w}],ListPlot[{{solx[[2,2,2]],solx[[1]]}},PlotStyle->{PointSize[0.03],Red}]}]

Fig .5 The second approach

Further steps lead to the following result, see Table 1.

Table 1.
The iteration steps of the McCormick envelopes approach

({{(Subscript[x, L]), (Subscript[x, U]), ->, (Subscript[x, solution]), w}, {-1, 2, ->, (1/2), -2}, {-1, (1/2), ->, -(1/4), -(1/2)}, {-(1/4), (1/2), ->, (1/8), -(1/8)}, {-(1/4), (1/8), ->, -(1/16), -(1/32)}, {-(1/16), (1/8), ->, (1/32), -(1/128)}, {-(1/16), (1/32), ->, -(1/64), -(1/512)}, {-(1/64), (1/32), ->, (1/128), -(1/2048)}})

1.3.2 Application of McCormic envelope

Let us employ McCormic envelopes approach to our simple quadratic problem. The new linear objective function of our example,

qL=27 (Subscript[w, 11])-18 (Subscript[w, 12]) +4 (Subscript[w, 22])-3 (Subscript[x, 2])

The box - type bounds are,

(Subscript[x, 1L])=-3;(Subscript[x, 1U])=4;(Subscript[x, 2L])=-2;(Subscript[x, 2U])=5;

The function values at the lower and upper bounds,

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->(Subscript[x, 1L]),(Subscript[x, 2])->(Subscript[x, 2L])}

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->(Subscript[x, 1U]),(Subscript[x, 2])->(Subscript[x, 2U])}

The additional inequality constraints, the McCormick envelopes are,

constraints=Apply[And,{(Subscript[w, 11])≥(Subscript[x, 1L]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1L])-(Subscript[x, 1L]) (Subscript[x, 1L]),(Subscript[w, 11])≥(Subscript[x, 1U]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1U])-(Subscript[x, 1U]) (Subscript[x, 1U]),
(Subscript[w, 11])≤(Subscript[x, 1U]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1L])-(Subscript[x, 1U]) (Subscript[x, 1L]),(Subscript[w, 11])≤(Subscript[x, 1])(Subscript[x, 1U])+(Subscript[x, 1L]) (Subscript[x, 1])-(Subscript[x, 1U]) (Subscript[x, 1L]),
(Subscript[w, 22])≥(Subscript[x, 2L]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2L])-(Subscript[x, 2L]) (Subscript[x, 2L]),(Subscript[w, 22])≥(Subscript[x, 2U]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2U])-(Subscript[x, 2U]) (Subscript[x, 2U]),
(Subscript[w, 22])≤(Subscript[x, 2U]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2L])-(Subscript[x, 2U]) (Subscript[x, 2L]),
(Subscript[w, 12])≥(Subscript[x, 1L]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2L])-(Subscript[x, 1L]) (Subscript[x, 2L]),(Subscript[w, 12])≥(Subscript[x, 1U]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2U])-(Subscript[x, 1U]) (Subscript[x, 2U]),
(Subscript[w, 12])≤(Subscript[x, 1U]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2L])-(Subscript[x, 1U]) (Subscript[x, 2L]),(Subscript[w, 12])≤(Subscript[x, 1])(Subscript[x, 2U])+(Subscript[x, 1L]) (Subscript[x, 2])-(Subscript[x, 2U]) (Subscript[x, 1L]),(Subscript[w, 11])>0,(Subscript[w, 12])≠0,(Subscript[w, 22])>0,
(Subscript[x, 1L])≤(Subscript[x, 1])≤(Subscript[x, 1U]),(Subscript[x, 2L])≤(Subscript[x, 2])≤(Subscript[x, 2U]),(Subscript[x, 1])≠0,(Subscript[x, 2])≠0,Element[{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])},Integers]}]

So we have a linear integer problem.

Contents

2.1 Nonlinear methods

2.1.1 Integer nonlinear optimization

2.1.2 Global numerical method

2.1.3 Local numerical method

2.2 Linear programming

2.3 Successive approximation

2.1 Nonlinear methods

2.1.1 Integer nonlinear optimization

To solve the linearized model, one may use nonlinear programming,

Minimize[{qL,constraints},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}]//AbsoluteTiming

where the minimum is

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->2,(Subscript[x, 2])->3}

Fig. 6. shows a new, reduced region of bounds using the solution as new upper limit.

Show[{p1,Graphics[{Yellow,Rectangle[{-3,-2},{4,5}]},Frame->True],Graphics[{Pink,Rectangle[{2,3},{-3,-2}]}],ListPlot[{{-3,-2},{4,5},{2,3}}->{"157","157","27"},LabelingFunction->Above,PlotStyle->PointSize[0.025]],p1}]

Fig .6 Solving linear problem, first approach

2.1.2 Global numerical method

Using numerical global minimization provides a considerable faster solution. Modifying input data structure,

p=Join[{Join[{qL},{constraints}]},{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}}]

then

NMinimize@@p//AbsoluteTiming

2.1.3 Local numerical method

Further improvement may be possible  via relaxing the non- equality constraints,

cons=constraints/.{(Subscript[x, 1])≠0->(Subscript[x, 1])>0,(Subscript[x, 2])≠0->(Subscript[x, 2])>0,(Subscript[w, 12])≠0->(Subscript[w, 12])>0};

and using local minimization,

p=Join[{Join[{qL},{cons}]},{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}}];

FindMinimum@@p//AbsoluteTiming

2.2 Linear programming

It goes without saying that the most proper method to solve the linear problem is the linear programming algorithm. We should construct the corresponding input.

The coefficient vector of the objective function,

c={0,-3,27,-18,4}

We introduce a small positive constant to exclude trivial solution,

ϵ=(10^(-3)); (*excluding the trivial solution {0,0}*)

Then the matrix form of the constraints,

m=(({{-2 (Subscript[x, 1L]) , 0, 1, 0, 0}, {-2(Subscript[x, 1U]), 0, 1, 0, 0}, {(Subscript[x, 1U])+(Subscript[x, 1L]), 0, -1, 0, 0}, {0, -2 (Subscript[x, 2L]) , 0, 0, 1}, {0, -2(Subscript[x, 2U]), 0, 0, 1}, {0, (Subscript[x, 2U])+(Subscript[x, 2L]), 0, 0, -1}, {-(Subscript[x, 2L]), -(Subscript[x, 1L]), 0, 1, 0}, {-(Subscript[x, 2U]), -(Subscript[x, 1U]), 0, 1, 0}, {(Subscript[x, 2L]), (Subscript[x, 1U]), 0, -1, 0}, {(Subscript[x, 2U]), (Subscript[x, 1L]), 0, -1, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 1}, {1, 0, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 1, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}}));b=(({{-(Subscript[x, 1L]) (Subscript[x, 1L])}, {-(Subscript[x, 1U]) (Subscript[x, 1U])}, {(Subscript[x, 1U]) (Subscript[x, 1L])}, {-(Subscript[x, 2L]) (Subscript[x, 2L])}, {-(Subscript[x, 2U]) (Subscript[x, 2U])}, {(Subscript[x, 2U]) (Subscript[x, 2L])}, {-(Subscript[x, 1L]) (Subscript[x, 2L])}, {-(Subscript[x, 1U]) (Subscript[x, 2U])}, {(Subscript[x, 1U]) (Subscript[x, 2L])}, {(Subscript[x, 2U]) (Subscript[x, 1L])}, {ϵ}, {ϵ}, {ϵ}, {(Subscript[x, 1L])}, {-(Subscript[x, 1U])}, {(Subscript[x, 2L])}, {-(Subscript[x, 2U])}, {ϵ}, {ϵ}}));

Now linear programming can be employed,

LinearProgramming[c,m,b//Flatten,Automatic,Integers]//Quiet//AbsoluteTiming

The running time for this toy problem can be longer than that of the local minimization, however, larger the problem, the linear programming is better.

2.3  Successive approximation

Using this  result {(Subscript[x, 1])=2, (Subscript[x, 2 ])= 3} as a new upper limit (see Fig.6), a new computation, as a second approach can be done,

(Subscript[x, 1L])=-3;(Subscript[x, 1U])=2;(Subscript[x, 2L])=-2;(Subscript[x, 2U])=3;

The values of the objective function at the lower bound,

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->(Subscript[x, 1L]),(Subscript[x, 2])->(Subscript[x, 2L])}

at the upper bound

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->(Subscript[x, 1U]),(Subscript[x, 2])->(Subscript[x, 2U])}

The new constraints are,

constraints=Apply[And,{(Subscript[w, 11])≥(Subscript[x, 1L]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1L])-(Subscript[x, 1L]) (Subscript[x, 1L]),(Subscript[w, 11])≥(Subscript[x, 1U]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1U])-(Subscript[x, 1U]) (Subscript[x, 1U]),
(Subscript[w, 11])≤(Subscript[x, 1U]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1L])-(Subscript[x, 1U]) (Subscript[x, 1L]),(Subscript[w, 11])≤(Subscript[x, 1])(Subscript[x, 1U])+(Subscript[x, 1L]) (Subscript[x, 1])-(Subscript[x, 1U]) (Subscript[x, 1L]),
(Subscript[w, 22])≥(Subscript[x, 2L]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2L])-(Subscript[x, 2L]) (Subscript[x, 2L]),(Subscript[w, 22])≥(Subscript[x, 2U]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2U])-(Subscript[x, 2U]) (Subscript[x, 2U]),
(Subscript[w, 22])≤(Subscript[x, 2U]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2L])-(Subscript[x, 2U]) (Subscript[x, 2L]),
(Subscript[w, 12])≥(Subscript[x, 1L]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2L])-(Subscript[x, 1L]) (Subscript[x, 2L]),(Subscript[w, 12])≥(Subscript[x, 1U]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2U])-(Subscript[x, 1U]) (Subscript[x, 2U]),
(Subscript[w, 12])≤(Subscript[x, 1U]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2L])-(Subscript[x, 1U]) (Subscript[x, 2L]),(Subscript[w, 12])≤(Subscript[x, 1])(Subscript[x, 2U])+(Subscript[x, 1L]) (Subscript[x, 2])-(Subscript[x, 2U]) (Subscript[x, 1L]),(Subscript[w, 11])>0,(Subscript[w, 12])≠0,(Subscript[w, 22])>0,
(Subscript[x, 1L])≤(Subscript[x, 1])≤(Subscript[x, 1U]),(Subscript[x, 2L])≤(Subscript[x, 2])≤(Subscript[x, 2U]),(Subscript[x, 1])≠0,(Subscript[x, 2])≠0,Element[{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])},Integers]}]

In array form,

m=(({{-2 (Subscript[x, 1L]) , 0, 1, 0, 0}, {-2(Subscript[x, 1U]), 0, 1, 0, 0}, {(Subscript[x, 1U])+(Subscript[x, 1L]), 0, -1, 0, 0}, {0, -2 (Subscript[x, 2L]) , 0, 0, 1}, {0, -2(Subscript[x, 2U]), 0, 0, 1}, {0, (Subscript[x, 2U])+(Subscript[x, 2L]), 0, 0, -1}, {-(Subscript[x, 2L]), -(Subscript[x, 1L]), 0, 1, 0}, {-(Subscript[x, 2U]), -(Subscript[x, 1U]), 0, 1, 0}, {(Subscript[x, 2L]), (Subscript[x, 1U]), 0, -1, 0}, {(Subscript[x, 2U]), (Subscript[x, 1L]), 0, -1, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 1}, {1, 0, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 1, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}}));b=(({{-(Subscript[x, 1L]) (Subscript[x, 1L])}, {-(Subscript[x, 1U]) (Subscript[x, 1U])}, {(Subscript[x, 1U]) (Subscript[x, 1L])}, {-(Subscript[x, 2L]) (Subscript[x, 2L])}, {-(Subscript[x, 2U]) (Subscript[x, 2U])}, {(Subscript[x, 2U]) (Subscript[x, 2L])}, {-(Subscript[x, 1L]) (Subscript[x, 2L])}, {-(Subscript[x, 1U]) (Subscript[x, 2U])}, {(Subscript[x, 1U]) (Subscript[x, 2L])}, {(Subscript[x, 2U]) (Subscript[x, 1L])}, {ϵ}, {ϵ}, {ϵ}, {(Subscript[x, 1L])}, {-(Subscript[x, 1U])}, {(Subscript[x, 2L])}, {-(Subscript[x, 2U])}, {ϵ}, {ϵ}}));

Then the result of the second approach is,

LinearProgramming[c,m,b//Flatten,Automatic,Integers]//Quiet//AbsoluteTiming

where the value of the objective function is,

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->1,(Subscript[x, 2])->2}

Then the region of the bounds can be further reduced, see Fig.7

Show[{p1,Graphics[{Yellow,Rectangle[{-3,-2},{4,5}]},Frame->True],Graphics[{Pink,Rectangle[{2,3},{-3,-2}]}],Graphics[{Green,Rectangle[{1,2},{-3,-2}]}],ListPlot[{{-3,-2},{4,5},{2,3},{1,2}}->{"157","157","27","    1"},LabelingFunction->Above,PlotStyle->PointSize[0.025]],p1}]

Fig .7 Solving linear problem, second approach

Alternatively, the result with global numerical minimization methods,

p=Join[{Join[{qL},{constraints}]},{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}}]

NMinimize@@p//AbsoluteTiming

or with local method

cons=constraints/.{(Subscript[x, 1])≠0->(Subscript[x, 1])>0,(Subscript[x, 2])≠0->(Subscript[x, 2])>0,(Subscript[w, 12])≠0->(Subscript[w, 12])>0};

p=Join[{Join[{qL},{cons}]},{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}}];

FindMinimum@@p//AbsoluteTiming

These last two method provided different solutions since {0, 1} ≠ {1, 2}. However the objective function has the same value in both solution points,

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->0,(Subscript[x, 2])->1}

We can continue the iteration with any of them, we shall get the same result.

Then the new lower limit for the next approach,

(Subscript[x, 1L])=1;(Subscript[x, 1U])=2;(Subscript[x, 2L])=2;(Subscript[x, 2U])=3;

The values of the objective function at the lower

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->(Subscript[x, 1L]),(Subscript[x, 2])->(Subscript[x, 2L])}

and the upper limit,

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->(Subscript[x, 1U]),(Subscript[x, 2])->(Subscript[x, 2U])}

see Fig.8.

The new constraints are,

constraints=Apply[And,{(Subscript[w, 11])≥(Subscript[x, 1L]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1L])-(Subscript[x, 1L]) (Subscript[x, 1L]),(Subscript[w, 11])≥(Subscript[x, 1U]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1U])-(Subscript[x, 1U]) (Subscript[x, 1U]),
(Subscript[w, 11])≤(Subscript[x, 1U]) (Subscript[x, 1])+(Subscript[x, 1])(Subscript[x, 1L])-(Subscript[x, 1U]) (Subscript[x, 1L]),(Subscript[w, 11])≤(Subscript[x, 1])(Subscript[x, 1U])+(Subscript[x, 1L]) (Subscript[x, 1])-(Subscript[x, 1U]) (Subscript[x, 1L]),
(Subscript[w, 22])≥(Subscript[x, 2L]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2L])-(Subscript[x, 2L]) (Subscript[x, 2L]),(Subscript[w, 22])≥(Subscript[x, 2U]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2U])-(Subscript[x, 2U]) (Subscript[x, 2U]),
(Subscript[w, 22])≤(Subscript[x, 2U]) (Subscript[x, 2])+(Subscript[x, 2])(Subscript[x, 2L])-(Subscript[x, 2U]) (Subscript[x, 2L]),
(Subscript[w, 12])≥(Subscript[x, 1L]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2L])-(Subscript[x, 1L]) (Subscript[x, 2L]),(Subscript[w, 12])≥(Subscript[x, 1U]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2U])-(Subscript[x, 1U]) (Subscript[x, 2U]),
(Subscript[w, 12])≤(Subscript[x, 1U]) (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 2L])-(Subscript[x, 1U]) (Subscript[x, 2L]),(Subscript[w, 12])≤(Subscript[x, 1])(Subscript[x, 2U])+(Subscript[x, 1L]) (Subscript[x, 2])-(Subscript[x, 2U]) (Subscript[x, 1L]),(Subscript[w, 11])>0,(Subscript[w, 12])≠0,(Subscript[w, 22])>0,
(Subscript[x, 1L])≤(Subscript[x, 1])≤(Subscript[x, 1U]),(Subscript[x, 2L])≤(Subscript[x, 2])≤(Subscript[x, 2U]),(Subscript[x, 1])≠0,(Subscript[x, 2])≠0,Element[{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])},Integers]}]

m=(({{-2 (Subscript[x, 1L]) , 0, 1, 0, 0}, {-2(Subscript[x, 1U]), 0, 1, 0, 0}, {(Subscript[x, 1U])+(Subscript[x, 1L]), 0, -1, 0, 0}, {0, -2 (Subscript[x, 2L]) , 0, 0, 1}, {0, -2(Subscript[x, 2U]), 0, 0, 1}, {0, (Subscript[x, 2U])+(Subscript[x, 2L]), 0, 0, -1}, {-(Subscript[x, 2L]), -(Subscript[x, 1L]), 0, 1, 0}, {-(Subscript[x, 2U]), -(Subscript[x, 1U]), 0, 1, 0}, {(Subscript[x, 2L]), (Subscript[x, 1U]), 0, -1, 0}, {(Subscript[x, 2U]), (Subscript[x, 1L]), 0, -1, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 0, 0, 1}, {1, 0, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}, {0, 1, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 1, 0, 0, 0}}));b=(({{-(Subscript[x, 1L]) (Subscript[x, 1L])}, {-(Subscript[x, 1U]) (Subscript[x, 1U])}, {(Subscript[x, 1U]) (Subscript[x, 1L])}, {-(Subscript[x, 2L]) (Subscript[x, 2L])}, {-(Subscript[x, 2U]) (Subscript[x, 2U])}, {(Subscript[x, 2U]) (Subscript[x, 2L])}, {-(Subscript[x, 1L]) (Subscript[x, 2L])}, {-(Subscript[x, 1U]) (Subscript[x, 2U])}, {(Subscript[x, 1U]) (Subscript[x, 2L])}, {(Subscript[x, 2U]) (Subscript[x, 1L])}, {ϵ}, {ϵ}, {ϵ}, {(Subscript[x, 1L])}, {-(Subscript[x, 1U])}, {(Subscript[x, 2L])}, {-(Subscript[x, 2U])}, {ϵ}, {ϵ}}));

and the new solution,

LinearProgramming[c,m,b//Flatten,Automatic,Integers]//Quiet//AbsoluteTiming

where the objective is

(qL/.{(Subscript[w, 11])->(Subscript[x, 1]) (Subscript[x, 1]),(Subscript[w, 12])->(Subscript[x, 1]) (Subscript[x, 2]),(Subscript[w, 22])->(Subscript[x, 2]) (Subscript[x, 2])})/.{(Subscript[x, 1]) ->1,(Subscript[x, 2])->3}

Show[{p1,Graphics[{Yellow,Rectangle[{-3,-2},{4,5}]},Frame->True],Graphics[{Pink,Rectangle[{2,3},{-3,-2}]}],Graphics[{Green,Rectangle[{1,2},{-3,-2}]}],Graphics[{Red,Rectangle[{1,2},{2,3}]}],ListPlot[{{-3,-2},{4,5},{2,3},{1,2},{1,3}}->{"157","157","27","    1","    0"},LabelingFunction->Above,PlotStyle->PointSize[0.025]],p1}]

Fig. 8 Solving linear problem, third approach

Employing alternative solution methods,

p=Join[{Join[{qL},{constraints}]},{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}}]

NMinimize@@p//AbsoluteTiming

cons=constraints/.{(Subscript[x, 1])≠0->(Subscript[x, 1])>0,(Subscript[x, 2])≠0->(Subscript[x, 2])>0,(Subscript[w, 12])≠0->(Subscript[w, 12])>0};

p=Join[{Join[{qL},{cons}]},{{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 22])}}];

FindMinimum@@p//AbsoluteTiming

No more approximation step is necessary, since the next iteration will give the same result, so it is a fix point of the iteration process.
The linear approach required roughly  ∼4 % of the nonlinear running time.

Contents

3.1 Converting mixed problem into a pure one

2.1  Converting mixed problem into a pure one

Until now, we have considered a pure integer problem. However, in case of mixed integer problem a part of the variables are continuous variables. This type of the problem can be transformed into a pure integer problem. Let us consider the following illustrative example.

Let the function to be maximized is,

u=3 (Subscript[x, 1])+5 (Subscript[x, 2]) +(Subscript[y, 1])+2 (Subscript[y, 2])

with continuous or “float” variables ((Subscript[y, 1]),(Subscript[y, 2])) ∈ ℛ, and with integer variables ((Subscript[x, 1]),(Subscript[x, 2])) ∈ 𝒵. 
First let us solve the continuous version of the problem. The constraints are

con=Apply[And,{2 (Subscript[x, 1])+(Subscript[x, 2])+(Subscript[y, 1])+2 (Subscript[y, 2])≤4,(Subscript[x, 1])+3 (Subscript[x, 2])+(Subscript[y, 1])+ (Subscript[y, 2])≤5,(Subscript[x, 1])≥0,(Subscript[x, 2])≥0,(Subscript[y, 1])≥0,(Subscript[y, 2])≥0}]

Then the continuous solution is

solC=NMaximize[{u,con},{(Subscript[x, 1]), (Subscript[x, 2]),(Subscript[y, 1]), (Subscript[y, 2])}]//Rationalize

Now, we introduce new integer variables as,

(Subscript[ξ, i])=(10^(Accuracy( (Subscript[w, i])))) (Subscript[y, i]),  i=1,2

In our case let Accuracy ((Subscript[w, i])) = 3

(* (Subscript[ξ, 1])=1000 (Subscript[y, 1]),(Subscript[ξ, 2])=1000 (Subscript[y, 2]) *)

In this way the continuous variable is considered as an integer one with 3 digit accuracy.

Then the objective with the new variables  is.

u=3 (Subscript[x, 1])+5 (Subscript[x, 2]) +(Subscript[ξ, 1])(10^(-3))+2 (Subscript[ξ, 2])(10^(-3))

where now all variables are integers,

con=Apply[And,{2 (Subscript[x, 1])+(Subscript[x, 2])+(Subscript[ξ, 1])(10^(-3))+2 (Subscript[ξ, 2])(10^(-3))≤4,(Subscript[x, 1])+3 (Subscript[x, 2])+(Subscript[ξ, 1])(10^(-3))+ (Subscript[ξ, 2])(10^(-3))≤5,(Subscript[x, 1])≥0,(Subscript[x, 2])≥0,(Subscript[ξ, 1])≥0,(Subscript[ξ, 2])≥0,Element[{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[ξ, 1]),(Subscript[ξ, 2])},Integers]}]

The solution is

solIC=Maximize[{u,con},{(Subscript[x, 1]), (Subscript[x, 2]),(Subscript[ξ, 1]), (Subscript[ξ, 2])}]//Rationalize

then

{(Subscript[y, 1]),(Subscript[y, 2])}=(10^(-3)){(Subscript[ξ, 1]),(Subscript[ξ, 2])}/.solIC[[2]]

(Subscript[y, 1])=.;(Subscript[y, 2])=.;

Contents

4.1 Problem description

4.2 Input data

4.3 Data preprocessing

4.4 Solution methods

4.4.1 Direct nonlinear programming

4.4.2 Continuous solution

4.4.3 A naive solution

4.4.4 Linearized problem

4.4.5 Iterative solution

4.4.6 Solution without iteration

4.4.7 Employing linear programming

4.1 Problem description

Highly accurate static GNSS positioning in surveying is achieved by the processing of relative phase ranges observed to the visible GNSS satellites at both the reference and the rover stations.To eliminate the time - dependent error sources such as the satellite and receiver clock error, the double differenced phase observations are formed and they are adjusted using a least squares adjustment. An alternative technique is to use extended Kalman - filter for this purpose.The observation equation of the double differenced phase observations has the following form:

(Subscript[ΔΔΦ, AB]^jk) =(Subscript[a, 1])(Subscript[δx, B]) + (Subscript[a, 2])(Subscript[δy, B])+(Subscript[a, 3])(Subscript[δz, B])+(Subscript[λδN, AB]^jk)

where (Subscript[ΔΔΦ, AB]^jk) is the double differences phase observations taken to the j-th and k-th satellite, (Subscript[δx, B]), (Subscript[δy, B]) and (Subscript[δz, B]) are the relative coordinate differences between the reference (A) and rover (B) stations, λ is the wavelength of the signal, (Subscript[δN, AB]^jk) is the double differenced phase ambiguity and j refers to the so-called pivot satellite, that is used as a reference for forming the double differences Let us assume that five satellites are measured concurrently on both the reference and the rover stations in two consecutive epochs. Since one satellite is used as a pivot satellite, four double differences are formed in each epoch. This means that altogether 8 observation equations are formed, which can be used to evaluate 7 unknowns (3 coordinate differences and 4 double-differenced phase ambiguities).
A usual solution of the problem is to estimate the unknowns using a least-squares adjustment, which provides a 'float' solution of the integer phase ambiguities. Consequently the computation of the integer least-squares estimates of the GPS cycle ambiguities reduces to a integer least-squares problem, see Teunissen (1995), Grafarend (2003) and Teunissen (2012).

(((y - A x - B z))^T)((Subscript[Q, y])^(-1))(y - A x - B z)->(Underscript[min
, x,z])

where y - vector of DD carrier phase observation in cycles,

A - design matrix for continuous valued parameters (baseline components,

B - design matrix for ambiguities,

x - unknown vector of continuous parameters, x ϵ ℛ

z - unknown ambiguity vector in cycles, z ϵ 𝒵

(Subscript[Q, y]) - weight matrix.

4.2 Input data

Let us assume that five satellites are measured concurrently on both the reference and the rover stations in two consecutive epochs. Since one satellite is used as a pivot satellite, four double differences are formed in each epoch. This means that altogether 8 observation equations are formed, which can be used to evaluate 7 unknowns (3 coordinate differences and 4 double-differenced phase ambiguities).

The actual values of the input arrays Rózsa (2017),

y=(({{2.883}, {3.020}, {2.586}, {4.405}, {2.922}, {3.121}, {2.663}, {4.267}}));A=(({{0.25, 0.20, 0.45}, {-0.30, 0.42, 0.56}, {-0.34, -0.20, 0.78}, {0.22, 0.54, 0.33}, {0.28, 0.23, 0.40}, {-0.31, 0.47, 0.60}, {-0.31, -0.22, 0.87}, {0.20, 0.48, 0.30}}));B=(({{0.19, 0, 0, 0}, {0, 0.19, 0, 0}, {0, 0, 0.19, 0}, {0, 0, 0, 0.19}, {0.19, 0, 0, 0}, {0, 0.19, 0, 0}, {0, 0, 0.19, 0}, {0, 0, 0, 0.19}}));

and the weight matrix is,

Q=IdentityMatrix[8];

4.3 Data preprocessing

First we rationalize the input data,

{yℛ,𝒬ℛ,𝒜ℛ,ℬℛ}=Rationalize[{y,Q,A,B}];

Now, in order to convert the problem pure integer programming, two new integer variables will be introduced to eliminate the “float” unknowns assuming three digits accuracy, see section 2.1,

(Subscript[ξ, 1])= 1000(Subscript[x, 1]),  (Subscript[ξ, 2])= 1000(Subscript[x, 2]),   (Subscript[ξ, 3])= 1000(Subscript[x, 3])

Then the objective function can be written as

objective=((yℛ-𝒜ℛ.{(Subscript[ξ, 1]) (10^(-3)),(Subscript[ξ, 2]) (10^(-3)),(Subscript[ξ, 3]) (10^(-3))}-ℬℛ.{(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]),(Subscript[z, 4])})//Flatten).Inverse[𝒬ℛ];

objective=objective.(yℛ-𝒜ℛ.{(Subscript[ξ, 1]) (10^(-3)),(Subscript[ξ, 2]) (10^(-3)),(Subscript[ξ, 3]) (10^(-3))}-ℬℛ.{(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]),(Subscript[z, 4])})//Simplify//First//Expand;

The number of the variables are

vars={(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]),(Subscript[z, 4]),(Subscript[ξ, 1]),(Subscript[ξ, 2]),(Subscript[ξ, 3])};

n=Length[vars];

4.4 Solution methods

4.4.1 Direct nonlinear programming

First, we try to solve the problem using brute force, via global nonlinear integer programming,

sol=Minimize[{objective,Element[vars,Integers]},vars]

It does not work. Let us employ the continuous solution.

4.4.2 Continuous solution

Global minimum on real field,

solR=NMinimize[objective,vars]

It is easy and fast.

4.4.3 A naive solution

Before applying the algorithm illustrated in the previous sections, a naive approach can be considered. Starting with constraints:  Floor[continuous result of variables] -1 ≤ variable≤  Round[real result of variables] + 1, the lower and upper limit for the variables are,

sB=Map[{Floor[#[[2]]]-1,Round[#[[2]]]+1}&,solR[[2]]]

then their constrains,

constraints=MapThread[#1[[1]]≤#2≤#[[2]]&,{sB,vars}]

Now the nonlinear integer minimization works,

Minimize[{objective,constraints,Element[vars,Integers]},vars]//AbsoluteTiming//N

Further improvement is possible by fixing the integer solutions and solving the reduced system for the continuous variables. Let us substitute the integer result into the objective,

solI=Take[%[[2,2]],{1,4}];

objectivE=objective/.solI;

Now we can solve the minimization problem only for the continuous variables using numerical global solver

solC=NMinimize[objectivE,{(Subscript[ξ, 1]),(Subscript[ξ, 2]),(Subscript[ξ, 3])}]

Then the solution of the mixed integer problem is,

solIC=Join[solI,solC[[2]]]

The value of the objective

objective/.solIC

This is a very simple relaxation technique but may unreliable and risky.

4.4.4  Linearized problem

Now, it is time to employ our algorithm demonstrated with the simple quadratic problem in the previous sections. First we compute the bounds. The matrix of the bilinear form is

Q=1/2 D[objective,{vars,2}];

The eigenvalues

λ=Eigenvalues[Q]//N

Rescaling with ((((10^(-3))))^2) since the introduction of the new integer variables, the ratio of the maximum  and minimum  values of λ,

κ=((Max[λ])/(Min[λ]))(10^(-6))

Then the radius of the hyper sphere with the real solution as center,

R=1/2 (Sqrt[n κ])

The equation of the sphere with the continuous solution as a center and with radius R,

eq=Total[Map[(((#[[1]]-#[[2]]))^2)&,solR[[2]]]]-(R^2)

Then the box - type bound is,

Map[Flatten[Solve[(eq/.Drop[solR[[2]],{#,#}])==0,vars[[#]]]]&,Range[n]]

and these coordinates will be the lower and upper bounds,

sB=Map[{Floor[#[[1]][[2]]],Round[#[[2]][[2]]]}&,%]

Now let us construct the input arrays for the linear programming.First we compute the terms which remains constant during the iteration procedure.

xLU=Table[{(Subscript[xL, i]),(Subscript[xU, i])},{i,1,n}]

The exponents of the variables in the objective function,

pow=Map[#[[1]]&,CoefficientRules[objective,vars]](*once*)

In order to employ McCormick envelope we select the quadratic terms resulting new variables (Subscript[w, i,j])

two=Select[pow,Max[#]==2&](*once*)

Now we compute the indexes of these variables. First,

Wii=Map[Flatten[{First[Flatten[Position[#,2]]],First[Flatten[Position[#,2]]]}]&,two](*once*)

and

oneone=Select[pow,And[Total[#]==2,Max[#]==1]&](*once*)

in addition

Wij=Map[Flatten[Position[#,1]]&,oneone](*once*)

then

W=Join[Wii,Wij](*once*)

The variables (Subscript[z, i]) having exponent one, remain

one=Select[pow,Total[#]==1&](*once*)

The part of the constraints independent on the variables of the objective

constraints3=Map[(Subscript[w, #[[1]],#[[2]]])>0&,W](*once*)

The name of the variables,

Vars=Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]](*once*)

The part of the objective function represented by the new variables,

objectivew=Map[Coefficient[objective,vars[[#[[1]]]] vars[[#[[2]]]]](Subscript[w, #[[1]],#[[2]]])&,W]

The other part of objective function with the original variables,

objectiveX=Table[First[Coefficient[objective,vars[[i]]]]vars[[i]],{i,1,n}]

The constant term in the objective can be left out,

term={First[CoefficientList[objective,vars[[1]]]]//First}

Then the new complete objective function,

objectivenew=Total[Join[objectivew,objectiveX,term]]

These terms are independent on the variables of the objective function and remain constant during the iteration.

4.4.5 Iterative solution

Now we can solve the linearized models via iteration. First we use nonlinear programming method.

1st iteration:

The actual lower and upper limits are,

xLU=sB;

Then constraints from the McCornick envelopes,

constraints1=Map[{(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[1]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[1]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[2]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[2]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[2]]]][[2]] vars[[#[[1]]]]+vars[[#[[2]]]]xLU[[#[[1]]]][[1]]-xLU[[#[[2]]]][[2]] xLU[[#[[1]]]][[1]]}&,W]//Flatten//Union;

The direct constraints for the original variables,

constraints2=MapThread[#2[[1]]≤#1≤#2[[2]]&,{vars,sB}]

Combining all constraints,

constraints=Apply[And,Join[constraints1,constraints2,constraints3,{Element[Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]],Integers]}]];

Solving via global nonlinear integer programming,

sol1=Minimize[{objectivenew,constraints},Vars]//AbsoluteTiming

The lower and upper bounds we used

sB

Changing the right hand side (upper limit), the new bounds,

sB={{-2,11},{0,12},{2,14},{4,16},{1085,1097},{1535,1546},{816,827}}

2nd iteration:

xLU=sB;

constraints1=Map[{(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[1]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[1]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[2]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[2]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[2]]]][[2]] vars[[#[[1]]]]+vars[[#[[2]]]]xLU[[#[[1]]]][[1]]-xLU[[#[[2]]]][[2]] xLU[[#[[1]]]][[1]]}&,W]//Flatten//Union;

constraints2=MapThread[#2[[1]]≤#1≤#2[[2]]&,{vars,sB}]

constraints=Apply[And,Join[constraints1,constraints2,constraints3,{Element[Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]],Integers]}]];

sol1=Minimize[{objectivenew,constraints},Vars]//AbsoluteTiming

The lower and upper bounds we used

sB

Now lower limit (left hand side) have been changed

sB={{6,11},{7,12},{8,14},{4,16},{1090,1097},{1540,1546},{822,827}}

3rd iteration:

xLU=sB;

constraints1=Map[{(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[1]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[1]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[2]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[2]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[2]]]][[2]] vars[[#[[1]]]]+vars[[#[[2]]]]xLU[[#[[1]]]][[1]]-xLU[[#[[2]]]][[2]] xLU[[#[[1]]]][[1]]}&,W]//Flatten//Union;

constraints2=MapThread[#2[[1]]≤#1≤#2[[2]]&,{vars,sB}]

constraints=Apply[And,Join[constraints1,constraints2,constraints3,{Element[Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]],Integers]}]];

sol1=Minimize[{objectivenew,constraints},Vars]//AbsoluteTiming

sB

sB={{9,11},{10,12},{11,14},{4,16},{1093,1097},{1543,1546},{825,827}}

4th iteration:

xLU=sB;

constraints1=Map[{(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[1]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[1]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[2]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[2]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[2]]]][[2]] vars[[#[[1]]]]+vars[[#[[2]]]]xLU[[#[[1]]]][[1]]-xLU[[#[[2]]]][[2]] xLU[[#[[1]]]][[1]]}&,W]//Flatten//Union;

constraints2=MapThread[#2[[1]]≤#1≤#2[[2]]&,{vars,sB}]

constraints=Apply[And,Join[constraints1,constraints2,constraints3,{Element[Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]],Integers]}]];

sol1=Minimize[{objectivenew,constraints},Vars]//AbsoluteTiming

sB

sB={{10,11},{11,12},{13,14},{4,16},{1095,1097},{1545,1546},{826,827}}

5th iteration:

xLU=sB;

constraints1=Map[{(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[1]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[1]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[2]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[2]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[2]]]][[2]] vars[[#[[1]]]]+vars[[#[[2]]]]xLU[[#[[1]]]][[1]]-xLU[[#[[2]]]][[2]] xLU[[#[[1]]]][[1]]}&,W]//Flatten//Union;

constraints2=MapThread[#2[[1]]≤#1≤#2[[2]]&,{vars,sB}]

constraints=Apply[And,Join[constraints1,constraints2,constraints3,{Element[Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]],Integers]}]];

sol1=Minimize[{objectivenew,constraints},Vars]//AbsoluteTiming

sB

sB={{10,11},{11,12},{13,14},{4,16},{1095,1097},{1545,1546},{826,827}}

Since the constraints did not change, we finish the iterations.

As post processing, fixing the integer solutions solving for the real variable

objectivE=objective/.{(Subscript[z, 1])->10,(Subscript[z, 2])->12,(Subscript[z, 3])->14,(Subscript[z, 4])->16}

NMinimize[objectivE,{(Subscript[ξ, 1]),(Subscript[ξ, 2]),(Subscript[ξ, 3])}]

Table 2

The steps of the iterations

({{iteration, (Subscript[z, 1L]), (Subscript[z, 1U]), (Subscript[z, 2L]), (Subscript[z, 2U]), (Subscript[z, 3L]), (Subscript[z, 3U]), (Subscript[z, 4L]), (Subscript[z, 4U]), (Subscript[ξ, 1L]), (Subscript[ξ, 1U]), (Subscript[ξ, 2L]), (Subscript[ξ, 2U]), (Subscript[ξ, 3L]), (Subscript[ξ, 3U])}, {0, -2, 21, 0, 23, 2, 25, 4, 27, 1085, 1108, 1535, 1558, 816, 840}, {1, -2, 11, 0, 12, 2, 14, 4, 16, 1085, 1097, 1535, 1546, 816, 827}, {2, 6, 11, 7, 12, 8, 14, 4, 16, 1090, 1097, 1540, 1546, 822, 827}, {3, 9, 11, 10, 12, 11, 14, 4, 16, 1093, 1097, 1543, 1546, 825, 827}, {4, 10, 11, 11, 12, 13, 14, 4, 16, 1095, 1097, 1545, 1546, 826, 827}, {5, 10, 11, 11, 12, 13, 14, 4, 16, 1095, 1097, 1545, 1546, 826, 827}})

ListPlot[{{11,6,9,10,10,10},{12,7,10,11,12,12},{14,8,11,13,14,14},{16,16,16,16,16,16}},Joined->True,PlotTheme->"Scientific",PlotLegends->{(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3]),(Subscript[z, 4])}]

Fig.9 The integer variables as function of iteration steps

ListPlot[{{97,90,93,95,97,97},{46,40,43,45,46,46},{27,22,25,26,26,26}},Joined->True,PlotTheme->"Scientific",PlotLegends->{(Subscript[ξ, 1])-1000,(Subscript[ξ, 2])-1500,(Subscript[ξ, 3])-800}]

Fig.10 The continuous  variables as function of iteration steps

4.4.6 Solution without iteration

One may reduce the computation time considerably to solve the linearized model directly without iteration with other methods than the symbolic nonlinear solver, Minimize. The starting data are

sB={{-2,21},{0,23},{2,25},{4,27},{1085,1108},{1535,1558},{816,840}}

xLU=sB;

constraints1=Map[{(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[1]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[1]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≥xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[2]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[2]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[1]]]][[2]] vars[[#[[2]]]]+vars[[#[[1]]]]xLU[[#[[2]]]][[1]]-xLU[[#[[1]]]][[2]] xLU[[#[[2]]]][[1]],(Subscript[w, #[[1]],#[[2]]])≤xLU[[#[[2]]]][[2]] vars[[#[[1]]]]+vars[[#[[2]]]]xLU[[#[[1]]]][[1]]-xLU[[#[[2]]]][[2]] xLU[[#[[1]]]][[1]]}&,W]//Flatten//Union;

constraints2=MapThread[#2[[1]]≤#1≤#2[[2]]&,{vars,sB}]

constraints=Apply[And,Join[constraints1,constraints2,constraints3,{Element[Join[vars,Map[(Subscript[w, #[[1]],#[[2]]])&,W]],Integers]}]];

p=Join[{Join[{objectivenew},{constraints}]},{Vars}];

Employing numerical global optimization,

NMinimize@@p//AbsoluteTiming

Alternatively numerical local optimization,

FindMinimum@@p//AbsoluteTiming

4.3.7 Employing linear programming

However the most adequate solution is the linear programming,

con1=Map[If[SameQ[#[[1]]≥#[[2]],#],#[[1]]-#[[2]]==0,#[[2]]-#[[1]]==0]&,constraints1];

con2=Map[{#[[3]]-#[[2]]==0,#[[2]]-#[[1]]==0}&,constraints2]//Flatten;

con3=Map[#[[1]]==0&,constraints3];

con=Join[con1,con2,con3];

c=Coefficient[objectivenew,Vars];

{b,A}=CoefficientArrays[con,Vars];

A1=Normal[A];b1=-Normal[b];

LinearProgramming[c,A1,b1,Automatic,Integers]//Quiet//AbsoluteTiming

We demonstrated that, the GPS phase ambiguity problem - a global quadratic, mixed integer programming problem - can be effectively solved via McCormick envelopes linearization and integer linear programming. Although, one should introduce new integer variables eliminating the continuous ones and compute the bounds for the global solution.

The author acknowledges the data of the GPS example for Dr.habil  Szabolcs Rózsa.

Teunissenn PJG, De Jonge PJ and Tiberius CCJM (1997).  The LAMBDA Method for Fast GPS Ambiguity Resolution. NAVIGATION Journal of the Institute of Navigation, Vol.44, No.3, p. 373 - 383, Wiley, USA

Teunissenn PJG, De Jonge PJ and Tiberius CCJM (1995).  The LAMBDA Method for Fast GPS Surveying, Presented at the International Symposium”GPS Technology Applications” Bucharest, Romania, September 26-29

https://www.researchgate.net/publication/246672882_THE_LAMBDA-METHOD_FOR_FAST_GPS_SURVEYING

Grafarend, E W (2003) Mixed Integer-Real Valued Adjustment (IRA) Problems: GPS Initial Cycle Ambiguity Resolution by Means of the LLL Algorithm, in Book 2003 Geodesy-The Challenge of the 3rd Millennium,
Editors: Professor Dr. Erik W. Grafarend, Dr. Friedrich W. Krumm, Dr. Volker S. Schwarze,
ISBN: 978-3-642-07733-3 (Print) 978-3-662-05296-9, 
(Online) https://link.springer.com/chapter/10.1007%2F978-3-662-05296-9_32

Teunissen P J G (2012) Towards a Unified Framework for GNSS Ambiguity Resolution: Problems &
Solutions, https://www.nav.ei.tum.de/fileadmin/w00bkq/www/Colloquium/colloquium_teunissen_slides.pdf

Rózsa Sz. (2017) Private Communication