diff --git a/src/math/simplex/model_based_opt.cpp b/src/math/simplex/model_based_opt.cpp
index abc93a0b9..e5db201a8 100644
--- a/src/math/simplex/model_based_opt.cpp
+++ b/src/math/simplex/model_based_opt.cpp
@@ -19,7 +19,7 @@ Revision History:
 --*/
 
 #include "model_based_opt.h"
-
+#include "uint_set.h"
 
 std::ostream& operator<<(std::ostream& out, opt::ineq_type ie) {
     switch (ie) {
@@ -94,7 +94,7 @@ namespace opt {
         unsigned_vector other;
         unsigned_vector bound_trail, bound_vars;
         while (!objective().m_vars.empty()) {
-            TRACE("opt", tout << "tableau\n";);
+            TRACE("opt", display(tout << "tableau\n"););
             var v = objective().m_vars.back();
             unsigned x = v.m_id;
             rational const& coeff = v.m_coeff;
@@ -111,7 +111,7 @@ namespace opt {
                 // => coeff*x <= -t2*coeff/a2
                 // objective + t2*coeff/a2 <= ub
 
-                mul_add(m_objective_id, - coeff/bound_coeff, bound_row_index);
+                mul_add(false, m_objective_id, - coeff/bound_coeff, bound_row_index);
                 m_rows[bound_row_index].m_alive = false;
                 bound_trail.push_back(bound_row_index);
                 bound_vars.push_back(x);
@@ -204,8 +204,13 @@ namespace opt {
         rational lub_val;
         rational const& x_val = m_var2value[x];
         unsigned_vector const& row_ids = m_var2row_ids[x];
+        uint_set visited;
         for (unsigned i = 0; i < row_ids.size(); ++i) {
             unsigned row_id = row_ids[i];
+            if (visited.contains(row_id)) {
+                continue;
+            }
+            visited.insert(row_id);
             row& r = m_rows[row_id];
             if (r.m_alive) {
                 rational a = get_coefficient(row_id, x);
@@ -219,13 +224,15 @@ namespace opt {
                         bound_row_index = row_id;
                         bound_coeff = a;
                     }
-                    else if ((is_pos && value < lub_val) || (!is_pos && value > lub_val)) {
+                    else if ((value == lub_val && r.m_type == opt::t_lt) ||
+                             (is_pos && value < lub_val) || 
+                             (!is_pos && value > lub_val)) {
                         other.push_back(bound_row_index);
                         lub_val = value;
                         bound_row_index = row_id;                            
                         bound_coeff = a;
                     }
-                    else if (bound_row_index != row_id) {
+                    else {
                         other.push_back(row_id);
                     }
                 }
@@ -253,11 +260,14 @@ namespace opt {
             }
             if (id < var_id) {
                 lo = mid + 1;
-            }
+            }			
             else {
-                hi = mid - 1;
+                hi = mid;
             }
         }
+        if (lo == r.m_vars.size()) {
+            return rational::zero();
+        }
         unsigned id = r.m_vars[lo].m_id;
         if (id == var_id) {
             return r.m_vars[lo].m_coeff;
@@ -293,17 +303,18 @@ namespace opt {
 
         SASSERT(a1 == get_coefficient(row_src, x));
         SASSERT(!a1.is_zero());
-   
+        SASSERT(row_src != row_dst);
+                
         if (m_rows[row_dst].m_alive) {
             rational a2 = get_coefficient(row_dst, x);
-            mul_add(row_dst, -a2/a1, row_src);            
+            mul_add(row_dst != m_objective_id && a1.is_pos() == a2.is_pos(), row_dst, -a2/a1, row_src);            
         }
     }
     
     //
     // set row1 <- row1 + c*row2
     //
-    void model_based_opt::mul_add(unsigned row_id1, rational const& c, unsigned row_id2) {
+    void model_based_opt::mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2) {
         if (c.is_zero()) {
             return;
         }
@@ -354,11 +365,12 @@ namespace opt {
         r1.m_coeff += c*r2.m_coeff;
         r1.m_vars.swap(m_new_vars);
         r1.m_value += c*r2.m_value;
-        if (r2.m_type == t_lt) {
+
+        if (!same_sign && r2.m_type == t_lt) {
             r1.m_type = t_lt;
         }
-        else if (r2.m_type == t_le && r1.m_type == t_eq) {
-            r1.m_type = t_le;
+        else if (same_sign && r1.m_type == t_lt && r2.m_type == t_lt) {
+            r1.m_type = t_le;        
         }
         SASSERT(invariant(row_id1, r1));
     }
diff --git a/src/math/simplex/model_based_opt.h b/src/math/simplex/model_based_opt.h
index 4efcf8f8f..b48a96edb 100644
--- a/src/math/simplex/model_based_opt.h
+++ b/src/math/simplex/model_based_opt.h
@@ -75,7 +75,7 @@ namespace opt {
 
         void resolve(unsigned row_src, rational const& a1, unsigned row_dst, unsigned x);
 
-        void mul_add(unsigned row_id1, rational const& c, unsigned row_id2);
+        void mul_add(bool same_sign, unsigned row_id1, rational const& c, unsigned row_id2);
 
         void set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, ineq_type rel);
         
diff --git a/src/test/model_based_opt.cpp b/src/test/model_based_opt.cpp
index 94ac81222..8c4148638 100644
--- a/src/test/model_based_opt.cpp
+++ b/src/test/model_based_opt.cpp
@@ -1,20 +1,121 @@
 #include "model_based_opt.h"
+#include "util.h"
+#include "uint_set.h"
 
 typedef opt::model_based_opt::var var;
 
-static void add_ineq(opt::model_based_opt& mbo, unsigned x, int a, unsigned y, int b, int k, opt::ineq_type rel) {
+static void add_ineq(opt::model_based_opt& mbo, unsigned x, int a, int k, opt::ineq_type rel) {
+    vector<var> vars;
+    vars.push_back(var(x, rational(a)));
+    mbo.add_constraint(vars, rational(k), rel);
+}
+
+static void add_ineq(opt::model_based_opt& mbo, 
+                     unsigned x, int a, 
+                     unsigned y, int b, int k, 
+                     opt::ineq_type rel) {
     vector<var> vars;
     vars.push_back(var(x, rational(a)));
     vars.push_back(var(y, rational(b)));
     mbo.add_constraint(vars, rational(k), rel);
 }
 
-static void add_ineq(opt::model_based_opt& mbo, unsigned x, int a, int k, opt::ineq_type rel) {
+static void add_ineq(opt::model_based_opt& mbo, 
+                     unsigned x, int a, 
+                     unsigned y, int b, 
+                     unsigned z, int c, int k, 
+                     opt::ineq_type rel) {
     vector<var> vars;
     vars.push_back(var(x, rational(a)));
+    vars.push_back(var(y, rational(b)));
+    vars.push_back(var(z, rational(c)));
     mbo.add_constraint(vars, rational(k), rel);
 }
 
+static void add_random_ineq(opt::model_based_opt& mbo,
+                            random_gen& r,
+                            svector<int>  const& values,
+                            unsigned max_vars,
+                            unsigned max_coeff) {
+    unsigned num_vars = values.size();
+    uint_set used_vars;
+    vector<var> vars;
+    int value = 0;
+    for (unsigned i = 0; i < max_vars; ++i) {
+        unsigned x = r(num_vars);
+        if (used_vars.contains(x)) {
+            continue;
+        }
+        used_vars.insert(x);
+        int coeff = r(max_coeff + 1);
+        if (coeff == 0) {
+            continue;
+        }
+        unsigned sign = r(2);
+		coeff = sign == 0 ? coeff : -coeff;
+        vars.push_back(var(x, rational(coeff)));
+        value += coeff*values[x];
+    }
+    unsigned abs_value = value < 0 ? - value : value;
+    // value + k <= 0
+    // k <= - value
+    // range for k is 2*|value|
+    // k <= - value - range
+    opt::ineq_type rel = opt::t_le;
+
+    int coeff = 0;
+    if (r(4) == 0) {
+        rel = opt::t_eq;
+        coeff = -value;
+    }
+    else {
+        if (abs_value > 0) {
+            coeff = -value - r(2*abs_value);
+        }
+        else {
+            coeff = 0;
+        }
+        if (coeff != -value && r(3) == 0) {
+            rel = opt::t_lt;
+        }   
+    }
+    mbo.add_constraint(vars, rational(coeff), rel);
+}
+                        
+static void check_random_ineqs(random_gen& r, unsigned num_vars, unsigned max_value, unsigned num_ineqs, unsigned max_vars, unsigned max_coeff) {
+    opt::model_based_opt mbo;
+    svector<int> values;
+    for (unsigned i = 0; i < num_vars; ++i) {
+        values.push_back(r(max_value + 1));
+        mbo.add_var(rational(values.back()));
+    }
+    for (unsigned i = 0; i < num_ineqs; ++i) {
+        add_random_ineq(mbo, r, values, max_vars, max_coeff);
+    }
+
+    vector<var> vars;
+    vars.reset();
+    vars.push_back(var(0, rational(2)));
+    vars.push_back(var(1, rational(-2)));
+    mbo.set_objective(vars, rational(0));
+
+    mbo.display(std::cout);
+    opt::inf_eps value = mbo.maximize();    
+    std::cout << "optimal: " << value << "\n";
+    mbo.display(std::cout);
+    for (unsigned i = 0; i < values.size(); ++i) {
+        std::cout << i << ": " << values[i] << " -> " << mbo.get_value(i) << "\n";
+    }
+}
+
+static void check_random_ineqs() {
+    random_gen r(1);
+
+    for (unsigned i = 0; i < 1009; ++i) {
+        check_random_ineqs(r, 4, 5, 5, 3, 6);
+    }
+}
+
 // test with upper bounds
 static void test1() {
     opt::model_based_opt mbo;
@@ -118,6 +219,7 @@ static void test4() {
 // test with mix of upper and lower bounds
 
 void tst_model_based_opt() {
+    check_random_ineqs();
     test1();
     test2();
     test3();