diff --git a/include/PetriEngine/PQL/Contexts.h b/include/PetriEngine/PQL/Contexts.h
index f9f56535e..c4d7e2e91 100644
--- a/include/PetriEngine/PQL/Contexts.h
+++ b/include/PetriEngine/PQL/Contexts.h
@@ -180,16 +180,18 @@ namespace PetriEngine {
 
             SimplificationContext(const MarkVal* marking,
                     const PetriNet* net, uint32_t queryTimeout, uint32_t lpTimeout,
-                    Simplification::LPCache* cache, uint32_t potencyTimeout = 0)
+                    Simplification::LPCache* cache, uint32_t potencyTimeout = 0, uint32_t printLevel = 0, bool useBigM = true)
                     : _queryTimeout(queryTimeout), _lpTimeout(lpTimeout),
                     _potencyTimeout(potencyTimeout) {
                 _negated = false;
+                _useBigM = useBigM;
                 _marking = marking;
                 _net = net;
                 _base_lp = buildBase();
                 _start = std::chrono::high_resolution_clock::now();
                 _cache = cache;
                 _markingOutOfBounds = false;
+                _printLevel = printLevel;
                 for(size_t i = 0; i < net->numberOfPlaces(); ++i) {
                     if (marking[i] >  std::numeric_limits<int32_t>::max()) { //too many tokens exceeding int32_t limits, LP solver will give wrong results
                         _markingOutOfBounds = true;
@@ -230,6 +232,14 @@ namespace PetriEngine {
 
             double getReductionTime();
 
+            uint32_t getPrintLevel() const{
+                return _printLevel;
+            }
+
+            bool useBigM() const {
+                return _useBigM;
+            }
+
             bool timeout() const {
                 auto end = std::chrono::high_resolution_clock::now();
                 auto diff = std::chrono::duration_cast<std::chrono::seconds>(end - _start);
@@ -254,10 +264,12 @@ namespace PetriEngine {
 
         private:
             bool _negated;
+            bool _useBigM;
             const MarkVal* _marking;
             bool _markingOutOfBounds;
             const PetriNet* _net;
             uint32_t _queryTimeout, _lpTimeout, _potencyTimeout;
+            uint32_t _printLevel;
             mutable glp_prob* _base_lp = nullptr;
             std::chrono::high_resolution_clock::time_point _start;
             Simplification::LPCache* _cache;
diff --git a/include/PetriEngine/Simplification/Member.h b/include/PetriEngine/Simplification/Member.h
index c30801b5b..a672d39b5 100644
--- a/include/PetriEngine/Simplification/Member.h
+++ b/include/PetriEngine/Simplification/Member.h
@@ -4,6 +4,7 @@
 #include <functional>
 #include <cassert>
 #include <cstring>
+#include <cstdint>
 //#include "../PQL/PQL.h"
 
 namespace PetriEngine {
diff --git a/include/PetriEngine/options.h b/include/PetriEngine/options.h
index 49f29a216..66a8c29c5 100644
--- a/include/PetriEngine/options.h
+++ b/include/PetriEngine/options.h
@@ -63,6 +63,8 @@ struct options_t {
     std::set<size_t> querynumbers;
     Strategy strategy = Strategy::DEFAULT;
     int queryReductionTimeout = 30, intervalTimeout = 10, partitionTimeout = 5, lpsolveTimeout = 10, initPotencyTimeout = 10;
+    bool lpUseBigM = true;
+    int lpPrintLevel = 0;
     TraceLevel trace = TraceLevel::None;
     bool use_query_reductions = true;
     uint32_t siphontrapTimeout = 0;
diff --git a/include/utils/structures/shared_string.h b/include/utils/structures/shared_string.h
index 3cc41398c..f2de0f2da 100644
--- a/include/utils/structures/shared_string.h
+++ b/include/utils/structures/shared_string.h
@@ -14,6 +14,7 @@
 #include <string>
 #include <memory>
 #include <algorithm>
+#include <cstdint>
 
 typedef const std::string const_string;
 typedef std::shared_ptr<const_string> shared_const_string;
diff --git a/src/PetriEngine/PQL/Contexts.cpp b/src/PetriEngine/PQL/Contexts.cpp
index f0fb0b404..a1f486878 100644
--- a/src/PetriEngine/PQL/Contexts.cpp
+++ b/src/PetriEngine/PQL/Contexts.cpp
@@ -153,6 +153,26 @@ namespace PetriEngine {
                     return nullptr;
                 }
             }
+
+            /*glp_add_rows(lp, _net->numberOfPlaces());
+            for (size_t p = 0; p < _net->numberOfPlaces(); p++) {
+                    std::vector<double> row(nCol+1, 0);
+                    std::vector<int> indices(1, 0);
+                    row.shrink_to_fit();
+                
+                    for (size_t t = 0; t < _net->numberOfTransitions(); t++) {
+                        if(_net->outArc(t, p) - _net->inArc(p, t) != 0){
+                            row[t+1]  = _net->outArc(t, p);
+                            row[t+1] -= _net->inArc(p, t);
+                            indices.push_back(t+1);
+                        }
+                    }
+
+                    glp_set_mat_row(lp, rowno, indices.size() - 1, indices.data(), row.data());
+                    glp_set_row_bnds(lp, rowno, GLP_LO, -1.0 * marking()[p], infty);
+                    ++rowno;
+            }*/
+            
             return lp;
         }
     }
diff --git a/src/PetriEngine/Simplification/LinearProgram.cpp b/src/PetriEngine/Simplification/LinearProgram.cpp
index 5733e549d..cd6330995 100644
--- a/src/PetriEngine/Simplification/LinearProgram.cpp
+++ b/src/PetriEngine/Simplification/LinearProgram.cpp
@@ -6,6 +6,7 @@
 #include "PetriEngine/Simplification/LinearProgram.h"
 #include "PetriEngine/Simplification/LPCache.h"
 #include "PetriEngine/PQL/Contexts.h"
+#include "VerifyPN.h"
 
 namespace PetriEngine {
     namespace Simplification {
@@ -46,6 +47,137 @@ namespace PetriEngine {
 
         constexpr auto infty = std::numeric_limits<REAL>::infinity();
 
+        std::string bound_or_inf(double bound){
+            std::string inf = "inf";
+            std::string ninf = "ninf";
+            return (std::fabs(bound) != infty)? std::to_string(bound) : ((bound > 0)? inf : ninf);
+        }
+
+        void printConstraints(const PQL::SimplificationContext& context, glp_prob* lp){
+            int nCols = glp_get_num_cols(lp);
+            int nRows = glp_get_num_rows(lp);
+            
+            std::vector<int> indices(nCols + 1);
+            std::vector<double> coef(nCols + 1);
+            int l;
+
+            for(int i = 1; i <= nRows; i++){
+                l = glp_get_mat_row(lp, i, indices.data(), coef.data());
+                std::cout << "Row " << i << "[len: " << l << "]: ";
+
+
+
+                bool first_print = true;
+                for(int j = 1; j <= l; j++){
+                    if(!first_print) std::cout << " + ";
+                    if((indices[j] - 1) < context.net()->numberOfTransitions()){
+                        std::cout << coef[j] << "*" << *context.net()->transitionNames()[indices[j] - 1];
+                    }else{
+                        std::cout << coef[j] << "*" << glp_get_col_name(lp, indices[j]);
+                    }
+                    first_print = false;
+                }
+
+                if(first_print){
+                    std::cout << "<empty-row>";
+                }
+                
+
+                double up = glp_get_row_ub(lp, i);
+                double lb = glp_get_row_lb(lp, i);
+                switch(glp_get_row_type(lp, i)){
+                    case GLP_FR:
+                        std::cout << " is unconstraint";
+                        break;
+                    case GLP_LO:
+                        std::cout << " >= " << bound_or_inf(lb);
+                        break;
+                    case GLP_UP:
+                        std::cout << " <= " << bound_or_inf(up);
+                        break;
+                    case GLP_DB:
+                        std::cout << " in [" << bound_or_inf(lb) << ", " << bound_or_inf(up) << "]";
+                        break;
+                    case GLP_FX:
+                        std::cout << " = " << bound_or_inf(lb);
+                        break;
+                }
+                std::cout << "\n";
+            }
+            
+        }
+        
+
+        void printSolution(const PQL::SimplificationContext& context, glp_prob* lp, bool is_mip){
+            if(context.getPrintLevel() == 0)
+                return;
+
+            if(context.getPrintLevel() > 1){
+                printConstraints(context, lp);
+            }
+            
+            int nCols = glp_get_num_cols(lp);
+            bool first_print = true;
+            if(is_mip){
+                std::cout << "MIP Solution: ";
+            }else{
+                std::cout << "LP Solution: ";
+            }
+            for(int i = 1; i <= nCols;i++){
+                double coef = (is_mip)? glp_mip_col_val(lp, i) : glp_get_col_prim(lp, i);
+                if(coef != 0){
+                    if(!first_print)
+                        std::cout << " + ";
+                    first_print = false;
+                    if( i - 1 < context.net()->numberOfTransitions()){
+                        std::cout << coef << "*" << *context.net()->transitionNames()[i - 1];
+                    }else{
+                        std::cout << coef << "*" << glp_get_col_name(lp, i);
+                    }
+                }
+            }
+
+            if(first_print)
+                std::cout << "<empty>\n";
+            else{
+                std::cout << "\n";
+            }
+        }
+
+        std::vector<Condition_ptr> bigMQuery(const PQL::SimplificationContext& context){
+            
+            std::vector<Condition_ptr> res;
+
+            for(int p = 0; p < context.net()->numberOfPlaces(); p++){
+                auto rhs = std::make_shared<IdentifierExpr>(context.net()->placeNames()[p]);
+                auto lhs = std::make_shared<LiteralExpr>(10000);
+                auto lt = std::make_shared<LessThanCondition>(lhs, rhs);
+                res.push_back(lt);
+            }
+
+            auto orQuery = std::make_shared<OrCondition>(res);
+            auto FQuery = std::make_shared<FCondition>(orQuery);
+            auto EFQuery = std::make_shared<ECondition>(FQuery);
+            std::vector<Condition_ptr> queries;
+            queries.push_back(EFQuery);
+            return queries;
+        }
+
+
+        bool checkBigMConstraint(const PQL::SimplificationContext& context){
+            return true;
+            auto queries = bigMQuery(context);
+            options_t options;
+            options.lpUseBigM = false;
+            simplify_queries(context.marking(), context.net(), queries, options, std::cout);
+            std::cout << "done big M\n";
+            if(queries[0]->isTriviallyFalse()){
+                return true;
+            }else{
+                return false;
+            }
+        }
+
         bool LinearProgram::isImpossible(const PQL::SimplificationContext& context, uint32_t solvetime) {
             bool use_ilp = true;
             auto net = context.net();
@@ -65,7 +197,7 @@ namespace PetriEngine {
             }
 
             const uint32_t nCol = net->numberOfTransitions();
-            const uint32_t nRow = net->numberOfPlaces() + _equations.size();
+            const uint32_t nRow = 2 * net->numberOfPlaces() + _equations.size();
 
             std::vector<REAL> row = std::vector<REAL>(nCol + 1);
             std::vector<int32_t> indir(std::max(nCol, nRow) + 1);
@@ -118,6 +250,108 @@ namespace PetriEngine {
                 glp_set_col_bnds(lp, i, GLP_LO, 0, infty);
             }
 
+            if(context.useBigM()){
+                for (size_t p = 0; p < net->numberOfPlaces(); p++) {      
+                        auto numColsLP = glp_get_num_cols(lp);
+                        std::vector<double> inRow(1, 0);
+                        std::vector<int> inIndices(1, 0);
+                        
+                        std::vector<double> outRow(1, 0);
+                        std::vector<int> outIndices(1, 0);
+                    
+
+                        bool has_loop = false;
+                        bool has_in   = false;
+                        //bool requires_tokens = false;
+
+                        //std::map<uint32_t, std::vector<uint32_t>> loops_by_weight;
+                        //int min_pos_loop = 0;
+                        uint32_t minimum_weight = UINT32_MAX;
+
+                        for (size_t t = 0; t < net->numberOfTransitions(); t++) {
+                            if(net->outArc(t, p) != 0 && net->inArc(p, t) != 0){
+                                has_loop = true;
+                                /*if(net->inArc(p, t) > context.marking()[p])
+                                    requires_tokens = true;
+                                
+                                if(net->outArc(t, p) > net->inArc(p, t) && net->inArc(p, t) < min_pos_loop)
+                                    min_pos_loop = net->inArc(p, t);*/
+
+                                minimum_weight = std::min(minimum_weight, net->inArc(p, t));
+                                
+                                outRow.push_back(1);
+                                outIndices.push_back(t+1);
+                            }
+                            else{
+                                if(net->inArc(p,t) != 0){
+                                    outRow.push_back(1);
+                                    outIndices.push_back(t+1);
+                                }
+                                if(net->outArc(t,p) != 0){
+                                    has_in = true;
+                                    
+                                    inRow.push_back(net->outArc(t, p));
+                                    inIndices.push_back(t+1);
+                                }
+                            }
+                        }
+
+                        if(!has_loop || !has_in || context.marking()[p] >= minimum_weight){
+                            //std::cout << "no loop applicable\n";
+                            continue;
+                        }
+                        else{
+                            //std::cout << "adding loop constraint for " << *net->placeNames()[p] <<"\n";
+                        }
+
+                        int idx = numColsLP + 1;
+                        glp_add_cols(lp, 2);
+
+                        for(int i = glp_get_num_cols(lp) - 1; i <= glp_get_num_cols(lp); i++){
+                            glp_set_obj_coef(lp, i, 1);
+                            glp_set_col_kind(lp, i, GLP_BV);
+                            //glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
+                        }
+
+                        const int OR_P = glp_get_num_cols(lp) ;
+                        std::string or_p_name = "NO_IMPLY_" + *net->placeNames()[p];
+                        glp_set_col_name(lp, OR_P, or_p_name.c_str());
+                        const int OR_Q = glp_get_num_cols(lp) - 1;
+                        std::string or_q_name = "YES_IMPLY_" + *net->placeNames()[p];
+                        glp_set_col_name(lp, OR_Q, or_q_name.c_str());
+                    
+                        const double needed_weight = (double) (minimum_weight - context.marking()[p]);
+                        inRow.push_back(needed_weight);
+                        inIndices.push_back(OR_P);
+
+                        outRow.push_back(-1.0 * 10000);
+                        outIndices.push_back(OR_Q);
+                    
+                        glp_add_rows(lp, 3);
+
+                        glp_set_mat_row(lp, rowno, inIndices.size() - 1, inIndices.data(), inRow.data());
+                        glp_set_row_bnds(lp, rowno, GLP_LO, needed_weight, infty);
+                        ++rowno;
+
+                        glp_set_mat_row(lp, rowno, outIndices.size() - 1, outIndices.data(), outRow.data());
+                        glp_set_row_bnds(lp, rowno, GLP_UP, 0.0, 0.0);
+                        ++rowno;
+
+                        double orRow[3];
+                        int orIndices[3];
+                        
+
+                        orRow[1] = 1.0;
+                        orRow[2] = 1.0;
+                        orIndices[1] = OR_P;
+                        orIndices[2] = OR_Q;
+
+                        glp_set_mat_row(lp, rowno, 2, orIndices, orRow);
+                        glp_set_row_bnds(lp, rowno, GLP_FX, 1.0, 1.0);
+                        ++rowno;
+                }
+            }
+           
             // Minimize the objective
             glp_set_obj_dir(lp, GLP_MIN);
             auto stime = glp_time();
@@ -138,36 +372,71 @@ namespace PetriEngine {
                 auto status = glp_get_status(lp);
                 if (status == GLP_OPT)
                 {
+                    
+                    #if 1
                     glp_iocp iset;
                     glp_init_iocp(&iset);
                     iset.msg_lev = 0;
+                    iset.tol_int = 1e-10;
                     iset.tm_lim = std::min<uint32_t>(std::max<uint32_t>(timeout - (stime - glp_time()), 1), 1000);
                     iset.presolve = GLP_OFF;
                     auto ires = glp_intopt(lp, &iset);
+                    #else
+                    glp_smcp settings_exact;
+                    settings_exact.msg_lev = 0;
+                    settings_exact.it_lim = INT_MAX;
+                    settings_exact.tm_lim = std::min<uint32_t>(std::max<uint32_t>(timeout - (stime - glp_time()), 1), 1000);
+                    auto ires = glp_exact(lp, &settings_exact);
+                    #endif
                     if (ires == GLP_ETMLIM)
                     {
+                        printSolution(context, lp, false);
                         _result = result_t::UKNOWN;
                         // std::cerr << "glpk mip: timeout" << std::endl;
                     }
                     else if (ires == 0)
-                    {
+                    {   
+                        #if 1
                         auto ist = glp_mip_status(lp);
+                        #else
+                        auto ist = glp_get_status(lp);
+                        #endif
                         if (ist == GLP_OPT || ist == GLP_FEAS || ist == GLP_UNBND) {
+                            printSolution(context, lp, true);
                             _result = result_t::POSSIBLE;
                         }
                         else
                         {
-                            _result = result_t::IMPOSSIBLE;
+                            if(context.useBigM()){
+                                if(checkBigMConstraint(context)){
+                                    _result = result_t::IMPOSSIBLE;
+                                }else{
+                                    _result = result_t::POSSIBLE;
+                                }
+                            }
+                            else{
+                                _result = result_t::IMPOSSIBLE;
+                            }
                         }
                     }
                 }
                 else if (status == GLP_FEAS || status == GLP_UNBND)
                 {
+                    printSolution(context, lp, false);
                     _result = result_t::POSSIBLE;
                 }
                 else
                 {
-                    _result = result_t::IMPOSSIBLE;
+                    if(context.useBigM()){
+                        if(checkBigMConstraint(context)){
+                            _result = result_t::IMPOSSIBLE;
+                        }else{
+                            _result = result_t::POSSIBLE;
+                        }
+                    }
+                    else{
+                        _result = result_t::IMPOSSIBLE;
+                    }
                 }
             }
             else if (result == GLP_ENOPFS || result == GLP_ENODFS || result == GLP_ENOFEAS)
@@ -363,7 +632,7 @@ namespace PetriEngine {
                 for (size_t i = 1; i <= nCol; i++) {
                     glp_set_obj_coef(tmp_lp, i, row[i]);
                     glp_set_col_kind(tmp_lp, i, GLP_IV);
-                    glp_set_col_bnds(tmp_lp, i, GLP_LO, 0, infty);
+                    glp_set_col_bnds(tmp_lp, i, GLP_DB, 0, (double) INT32_MAX);
                 }
 
                 auto rs = glp_simplex(tmp_lp, &settings);
diff --git a/src/PetriEngine/options.cpp b/src/PetriEngine/options.cpp
index f5379408d..7bc04f96e 100644
--- a/src/PetriEngine/options.cpp
+++ b/src/PetriEngine/options.cpp
@@ -346,6 +346,13 @@ bool options_t::parse(int argc, const char** argv) {
             if (sscanf(argv[++i], "%d", &queryReductionTimeout) != 1 || queryReductionTimeout < 0) {
                 throw base_error("Argument Error: Invalid query reduction timeout argument ", std::quoted(argv[i]));
             }
+        }else if (std::strcmp(argv[i], "-lpp") == 0 || std::strcmp(argv[i], "--lp-print") == 0) {
+            if (i == argc - 1) {
+                throw base_error("Missing number after ", std::quoted(argv[i]));
+            }
+            if (sscanf(argv[++i], "%d", &lpPrintLevel) != 1 || lpPrintLevel < 0) {
+                throw base_error("Argument Error: Invalid lp print level argument ", std::quoted(argv[i]));
+            }
         } else if (std::strcmp(argv[i], "--init-potency-timeout") == 0) {
             if (i == argc - 1) {
                 throw base_error("Missing number after ", std::quoted(argv[i]));
diff --git a/src/VerifyPN.cpp b/src/VerifyPN.cpp
index c669c2990..3e7c0513e 100644
--- a/src/VerifyPN.cpp
+++ b/src/VerifyPN.cpp
@@ -388,6 +388,7 @@ std::vector<Condition_ptr> getLTLQueries(const std::vector<Condition_ptr>& ctlSt
     return ltlQueries;
 }
 
+
 #ifdef VERIFYPN_MC_Simplification
 std::mutex spot_mutex;
 #endif
@@ -578,7 +579,7 @@ void simplify_queries(const MarkVal* marking,
                     if (options.logic == TemporalLogic::LTL) {
                         if (options.queryReductionTimeout == 0 || qt == 0) continue;
                         SimplificationContext simplificationContext(marking, net, qt,
-                            options.lpsolveTimeout, &cache);
+                            options.lpsolveTimeout, &cache, 0, options.lpPrintLevel, options.lpUseBigM);
                         if (simplificationContext.markingOutOfBounds()) {
                             std::cout << "WARNING: Initial marking contains a place or places with too many tokens. Query simplifaction for LTL is skipped.\n";
                             break;
@@ -609,7 +610,7 @@ void simplify_queries(const MarkVal* marking,
 
                     if (options.queryReductionTimeout > 0 && qt > 0) {  
                         SimplificationContext simplificationContext(marking, net, qt,
-                            options.lpsolveTimeout, &cache);
+                            options.lpsolveTimeout, &cache, 0, options.lpPrintLevel, options.lpUseBigM);
                         if (simplificationContext.markingOutOfBounds()) {
                             std::cout << "WARNING: Initial marking contains a place or places with too many tokens. Query simplifaction is skipped.\n";
                             break;