Added setup.py file. Added usage description README

Author: bstellato

.gitignore

@@ -99,3 +99,6 @@ ENV/
 
 # Backups
 results_bak/
+
+# Vim backup files
+*.swp

README.md

@@ -1,7 +1,60 @@
-# MIQP python solver based on OSQP
-The solver lives in the `miosqp` package. In order to run it you need to install the [OSQP Python interface](https://github.com/bstellato/osqp).
+# Mixed-Integer Quadratic Program Python Solver Based on OSQP
 
-## Run example
+miOSQP solves an mixed-integer quadratic programs (MIQPs) of the form
+```
+minimize        0.5 x' P x + q' x
+
+subject to      l <= A x <= u
+                x[i] in Z for i in i_idx
+                i_l[i] <= x[i] <= i_u[i] for i in i_idx
+```
+where `i_idx` is a vector of indices of which variables are integer and `i_l`, `i_u` are the lower and upper bounds on the integer variables respectively.
+
+
+## Installation
+To install the package simply run
+```
+python setup.py install
+```
+it depends on [OSQP](osqp.readthedocs.io), numpy and scipy.
+
+
+## Usage
+To solve a MIQP we need to run
+```python
+import miosqp
+m = miosqp.MIOSQP()
+m.setup(P, q, A, l, u, i_idx, i_l, i_u)
+results = m.solve()
+```
+where `P` is a symmetric positive semidefinite matrix and `A` a matrix. 
+`P` and `A` are both in the scipy sparse CSC format.
+
+The returned object `results` contains 
+-  `x`: the solution
+-  `upper_glob`: the cost function upper bound
+-  `run_time`: the solution time
+-  `status`: the status
+-  `osqp_solve_time`: the OSQP solve time as a percentage of the total solution time
+-  `osqp_iter_avg`: the OSQP average number of iterations for each QP sub-problem solution
+
+
+### Update problem vectors
+Problem vectors can be updated without running the setup again. It can be done with
+```python
+m.update_vectors(q=q_new, l=l_new, u=u_new)
+```
+
+### Set initial solution guess
+The initial guess can speedup the branch-and-bound algorithm significantly.
+To set an initial feasible solution `x0` we can run
+```python
+m.set_x0(x0)
+```
+
+
+
+## Run examples
 In order to run the examples from to compare with GUROBI, after installing the python insterface, you need to install [mathprogbasepy](https://github.com/bstellato/mathprogbasepy). Examples can be found in the `examples` folder.
 
 -   Random MIQPs
@@ -10,5 +63,5 @@ In order to run the examples from to compare with GUROBI, after installing the p
 Note that you need [pandas](http://pandas.pydata.org/) package for storing the results dataframe and [tqdm](https://github.com/noamraph/tqdm) package for the progress bar.
 
 
-## Maximum number of iterations
-For some problem instances, OSQP reaches the maximum number of iterations. In order to deal with them, they are dumped to different files in the `max_iter_examples` folder. In order to load them separately and solve them with OSQP, you can run `examples/run_maxiter_problem.py` file.
+<!-- ## Maximum number of iterations -->
+<!-- For some problem instances, OSQP reaches the maximum number of iterations. In order to deal with them, they are dumped to different files in the `max_iter_examples` folder. In order to load them separately and solve them with OSQP, you can run `examples/run_maxiter_problem.py` file. -->

examples/power_converter/power_converter.py

@@ -462,7 +462,7 @@ def compute_mpc_input(self, x0, u_prev, solver='gurobi'):
                                  'eps_rel': 1e-03,
                                  'eps_prim_inf': 1e-04,
                                  #  'rho': 0.001,
-                                 'rho': 0.1,
+                                 #  'rho': 0.1,
                                  'verbose': False}
                 self.solver = miosqp.MIOSQP()
                 self.solver.setup(qp.P, q, qp.A, qp.l,
@@ -490,21 +490,22 @@ def compute_mpc_input(self, x0, u_prev, solver='gurobi'):
                 # import ipdb; ipdb.set_trace()
 
 
-
             if res_miosqp.status != miosqp.MI_SOLVED:
                 import ipdb; ipdb.set_trace()
             u = res_miosqp.x
             obj_val = res_miosqp.upper_glob
             solve_time = res_miosqp.run_time
-
+            osqp_solve_time = 100 * res_miosqp.osqp_solve_time / res_miosqp.run_time
 
         # Get first input
         u0 = u[:6]
 
         if solver == 'miosqp':
-            return u0, obj_val, solve_time, u, res_miosqp.osqp_iter_avg
+            return u0, obj_val, solve_time, u, \
+                    osqp_solve_time, \
+                    res_miosqp.osqp_iter_avg
         else:
-            return u0, obj_val, solve_time, u, 0
+            return u0, obj_val, solve_time, u, 0, 0
 
     def simulate_one_step(self, x, u):
         """
@@ -515,7 +516,6 @@ def simulate_one_step(self, x, u):
 
         return xnew, ynew
 
-
     def compute_signals(self, X):
         """
         Compute signals for plotting
@@ -533,7 +533,6 @@ def compute_signals(self, X):
         for i in range(T_final):
             Y_star_phase[:, i] = self.params.invP.dot(X[4:6, i])
 
-
         # Compute torque
         T_e = np.zeros(T_final)
         for i in range(T_final):
@@ -601,6 +600,7 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
         if solver == 'miosqp':
             # If miosqp, set avg numer of iterations to 0
             miosqp_avg_osqp_iter = 0
+            miosqp_osqp_avg_time = 0
 
         # Rename some variables for notation ease
         nx = self.dyn_system.A.shape[0]
@@ -629,7 +629,7 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
         for i in tqdm(range(T_final)):
 
             # Compute mpc inputs
-            U[:, i], obj_vals[i], time_temp, u_prev, osqp_iter = \
+            U[:, i], obj_vals[i], time_temp, u_prev, osqp_time, osqp_iter = \
                 self.compute_mpc_input(X[:, i], u_prev, solver=solver)
 
             # Store time if after the init periods
@@ -646,10 +646,13 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
             if solver == 'miosqp':
                 # Append average number of osqp iterations
                 miosqp_avg_osqp_iter += osqp_iter
+                miosqp_osqp_avg_time += osqp_time
 
         if solver == 'miosqp':
-            # Divide total number of average iterations by time steps
+            # Divide total number of average OSQP iterations 
+            # and solve time by time steps
             miosqp_avg_osqp_iter /= T_final
+            miosqp_osqp_avg_time /= T_final
 
         # Compute additional signals for plotting
         Y_phase, Y_star_phase, T_e, T_e_des = self.compute_signals(X)
@@ -667,5 +670,6 @@ def simulate_cl(self, N, steady_trans, solver='gurobi', plot=False):
 
         if solver == 'miosqp':
             stats.miosqp_avg_osqp_iter = miosqp_avg_osqp_iter
+            stats.miosqp_osqp_avg_time = miosqp_osqp_avg_time
 
         return stats

examples/power_converter/run_example.py

@@ -86,73 +86,60 @@ def run_example():
     miosqp_min_time = np.zeros(len(N_adp))
     miosqp_max_time = np.zeros(len(N_adp))
     miosqp_avg_osqp_iter = np.zeros(len(N_adp))
-
+    miosqp_osqp_avg_time = np.zeros(len(N_adp))
 
     # Simulate model
-
     for i in range(len(N_adp)):
 
-        stats_gurobi = model.simulate_cl(N_adp[i], flag_steady_trans, solver='gurobi')
+        stats_gurobi = model.simulate_cl(N_adp[i], flag_steady_trans, 
+                                         solver='gurobi')
 
         gurobi_std_time[i] = stats_gurobi.std_solve_time
         gurobi_avg_time[i] = stats_gurobi.avg_solve_time
         gurobi_min_time[i] = stats_gurobi.min_solve_time
         gurobi_max_time[i] = stats_gurobi.max_solve_time
 
-
         # Make plots for horizon 3
         if N_adp[i] == 3:
             plot_flag = 1
         else:
             plot_flag = 0
 
-        stats_miosqp = model.simulate_cl(N_adp[i], flag_steady_trans, solver='miosqp', plot=plot_flag)
-
+        stats_miosqp = model.simulate_cl(N_adp[i], flag_steady_trans, 
+                                         solver='miosqp', plot=plot_flag)
 
         miosqp_std_time[i] = stats_miosqp.std_solve_time
         miosqp_avg_time[i] = stats_miosqp.avg_solve_time
         miosqp_min_time[i] = stats_miosqp.min_solve_time
         miosqp_max_time[i] = stats_miosqp.max_solve_time
         miosqp_avg_osqp_iter[i] = stats_miosqp.miosqp_avg_osqp_iter
-
-
+        miosqp_osqp_avg_time[i] = stats_miosqp.miosqp_osqp_avg_time
 
     # Create pandas dataframe
-    timings = pd.DataFrame({ 'grb_avg' : gurobi_avg_time,
-                             'grb_std' : gurobi_std_time,
-                             'grb_min' : gurobi_min_time,
-                             'grb_max' : gurobi_max_time,
-                             'miosqp_avg' : miosqp_avg_time,
-                             'miosqp_std' : miosqp_std_time,
-                             'miosqp_min' : miosqp_min_time,
-                             'miosqp_max' : miosqp_max_time,
-                             'miosqp_avg_osqp_iter': miosqp_avg_osqp_iter},
-                             index=N_adp)
-
-
-
+    timings = pd.DataFrame({'T': N_adp,
+                            'grb_avg': gurobi_avg_time,
+                            'grb_std': gurobi_std_time,
+                            'grb_min': gurobi_min_time,
+                            'grb_max': gurobi_max_time,
+                            'miosqp_avg': miosqp_avg_time,
+                            'miosqp_std': miosqp_std_time,
+                            'miosqp_min': miosqp_min_time,
+                            'miosqp_max': miosqp_max_time,
+                            'miosqp_osqp_avg_time': miosqp_osqp_avg_time,
+                            'miosqp_avg_osqp_iter': miosqp_avg_osqp_iter})
 
     print("Results")
     print(timings)
 
     timings.to_csv('results/power_converter_timings.csv')
 
-    # # Converting results to latex table and storing them to a file
-    # formatter = lambda x: '%1.2f' % x
-    # latex_table = timings.to_latex(header=False,
-    #                                float_format=formatter)
-    # f = open('results/power_converter_miqp.tex', 'w')
-    # f.write(latex_table)
-    # f.close()
-
-
     # Create error plots with fill_between
     plt.figure()
     ax = plt.gca()
     plt.semilogy(N_adp, gurobi_avg_time, color=colors['o'],
-             label='GUROBI')
+                 label='GUROBI')
     plt.semilogy(N_adp, miosqp_avg_time, color=colors['b'],
-             label='miOSQP')
+                 label='miOSQP')
     plt.xticks(N_adp)
     ax.set_xlabel(r'Horizon length $T$')
     ax.set_ylabel(r'Time [s]')
@@ -162,8 +149,3 @@ def run_example():
     # plt.show()
     plt.savefig('results/power_converter_timings.pdf')
 
-
-
-
-
-    #