Ill-Conditioned Constraints and Solver Performance Issues

Hi all,

I'm working on an optimization model using Gurobi (via JuMP in Julia), and I'm experiencing solver performance that is highly sensitive to small changes in input data. In particular, solve times are much longer than expected and can vary dramatically with minor perturbations in the data.

Here is the log:

Optimize a model with 616688 rows, 14687 columns and 455711 nonzeros
Model fingerprint: 0x8092bf4e
Variable types: 14232 continuous, 455 integer (455 binary)
Coefficient statistics:
 Matrix range     [1e-13, 3e+06]
 Objective range  [1e-02, 1e+00]
 Bounds range     [2e-02, 3e-02]
 RHS range        [3e-02, 4e+02]
Warning: Model contains large matrix coefficient range

User MIP start did not produce a new incumbent solution

Presolve removed 586971 rows and 2395 columns
Presolve time: 0.06s
Presolved: 29717 rows, 12292 columns, 83227 nonzeros
Variable types: 12270 continuous, 22 integer (22 binary)
Performing another presolve...
Presolve removed 18851 rows and 7614 columns
Presolve time: 0.12s
Found heuristic solution: objective 0.0257797
Root relaxation presolve removed 18 rows and 18 columns
Root relaxation presolved: 10848 rows, 4660 columns, 32860 nonzeros

Extra simplex iterations after uncrush: 7

Root relaxation: objective 1.651929e-02, 8899 iterations, 1.12 seconds (2.01 work units)

   Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

    0     0    0.01652    0   14    0.02578    0.01652  35.9%     -    2s
    0     0    0.01652    0   14    0.02578    0.01652  35.9%     -    2s
H    0     0                       0.0255183    0.01857  27.2%     -    2s
    0     2    0.01857    0   13    0.02552    0.01857  27.2%     -    3s
    3     6    0.02300    2   12    0.02552    0.02089  18.2% 13871   25s
    7     7    0.02140    3   15    0.02552    0.02110  17.3% 21771   32s
*    9     7               3       0.0254660    0.02110  17.1% 17611   32s
   12     8    0.02357    4    2    0.02547    0.02113  17.0% 13292   40s
H   15    12                       0.0254316    0.02113  16.9% 10706   63s
*   20    12               5       0.0254259    0.02209  13.1% 23823   63s
   23    15    0.02414    5    2    0.02543    0.02215  12.9% 20753   95s
   30    16    0.02449    6   12    0.02543    0.02243  11.8% 23252  129s
H   32    15                       0.0248584    0.02269  8.74% 26956  129s
   35    21    0.02393    6   12    0.02486    0.02270  8.67% 26565  156s
   43    22    0.02273    8   11    0.02486    0.02273  8.55% 26841  171s
*   46    22              10       0.0243274    0.02273  6.56% 25705  171s
*   50    22               9       0.0235045    0.02273  3.29% 23659  171s
   54     8  postponed   10         0.02350    0.02273  3.29% 23149  185s
   64     9     cutoff    8         0.02350    0.02275  3.23% 20101  212s
   77     4     cutoff   15         0.02350    0.02278  3.06% 18102  221s

Explored 100 nodes (1488174 simplex iterations) in 221.87 seconds (447.89 work units)
Thread count was 8 (of 20 available processors)

Solution count 8: 0.0235045 0.0243274 0.0248584 ... 0.0257797

Optimal solution found (tolerance 0.00e+00)
Warning: max constraint violation (2.8376e-09) exceeds tolerance
Best objective 2.350452626891e-02, best bound 2.350452626891e-02, gap 0.0000%

My model includes two key constraints of the form:
for k in 1:λ
   @constraint(model, real(Y) * psi_r[:,k] - imag(Y) * psi_i[:,k] .== A * real(I[:,k]))
   @constraint(model, real(Y) * psi_i[:,k] + imag(Y) * psi_r[:,k] .== A * imag(I[:,k]))
end

Both Y and I are fixed input data.  A is a binary variable matrix. psi_r and psi_i are continuous variables. The issue is that:

• The nonzero entries in |I| range from 1e-13 to 1.
• The nonzero entries in |Y| range from 0.01 to 3e6.
• Y  have a very large condition number.

I suspect that this poor numerical conditioning is the root of the performance problem. 
I’m trying to understand

1. Do the attached solver logs suggest any additional problems besides numerical instability?
2. What is the best way to handle numerical instability in my constraints?

Any help would be greatly appreciated. Thanks!