don't warmstart the lpec with previous solution

examples/cardinality_optimization.m

@@ -9,8 +9,8 @@
 import casadi.*
 
 %% Problem formulation
-n = 300; % number of variabples
-m = 150;
+n = 25; % number of variabples
+m = 10;
 A = rand(m,n);
 b = rand(m,1);
 c = 1-2*rand(n,1);
@@ -43,32 +43,44 @@
 solver_initalization = struct('x0', x0,'lbx',lbx,'ubx',ubx,'lbg',lbg,'ubg',ubg);
 settings = HomotopySolverOptions();
 settings.homotopy_parameter_steering = "Direct";
-[result_homotopy,stats_homotopy] = mpec_homotopy_solver(mpec,solver_initalization,settings);
-f_opt_homotopy = full(result_homotopy.f);
-w_opt_homotopy = full(result_homotopy.x);
-%%  Settings
+[result_reg,stats_homotopy] = mpec_homotopy_solver(mpec,solver_initalization,settings);
+f_opt_reg = full(result_reg.f);
+w_opt_reg = full(result_reg.x);
+%%  mpecopt settings
 solver_settings = mpecopt.Options();
 % solver_settings.settings_lpec.lpec_solver = "Highs_casadi";
-% solver_settings.settings_lpec.lpec_solver = "Gurobi";
-solver_settings.rho_TR_phase_i_init = 1e-1;
-solver_settings.relax_and_project_homotopy_parameter_steering = "Direct";
-% solver_settings.initalization_strategy = "TakeInitialGuessDirectly";
+solver_settings.settings_lpec.lpec_solver = "Gurobi";
+% solver_settings.rho_TR_phase_i_init = 1e-1;
+% solver_settings.rho_TR_phase_ii_init = 1e-1;
+% solver_settings.relax_and_project_homotopy_parameter_steering = "Direct";
+% solver_settings.initialization_strategy = "FeasibilityEllInfGeneral";
 
+solver_initalization = struct('x0', x0, 'lbx',lbx, 'ubx',ubx,'lbg',lbg, 'ubg',ubg);
+solver = mpecopt.Solver(mpec, solver_settings);
+[result_active_set,stats_active_set] = solver.solve(solver_initalization);
+
+%
 solver_initalization = struct('x0', x0, 'lbx',lbx, 'ubx',ubx,'lbg',lbg, 'ubg',ubg);
 [result_active_set,stats_active_set] = mpec_optimizer(mpec, solver_initalization, solver_settings);
+% solver = mpecopt.Solver(mpec, solver_settings);
+
 %%
 w_opt_active_set = full(result_active_set.x);
 f_opt_active_set = full(result_active_set.f);
+
+x_opt_active_set = w_opt_active_set(1:n);
+x_opt_reg = w_opt_reg(1:n);
+cardinality_active_set = sum(heaviside(abs(x_opt_active_set)-1e-3));
+cardinality_reg  = sum(heaviside(abs(x_opt_reg)-1e-3));
+
+fprintf('\n Cardinality reg: %d, Cardinality mpecopt: %d \n',cardinality_reg,cardinality_active_set)
 fprintf('\n-------------------------------------------------------------------------------\n');
 fprintf('Method \t\t Objective \t comp_res \t n_biactive \t CPU time (s)\t Sucess\t Stat. type\n')
 fprintf('-------------------------------------------------------------------------------\n');
-fprintf('homotopy \t %2.2e \t %2.2e \t\t %d \t\t\t %2.2f \t\t\t\t %d\t %s\n',f_opt_homotopy,stats_homotopy.comp_res,stats_homotopy.n_biactive,stats_homotopy.cpu_time_total,stats_homotopy.success,stats_homotopy.multiplier_based_stationarity)
+fprintf('Reg \t %2.2e \t %2.2e \t\t %d \t\t\t %2.2f \t\t\t\t %d\t %s\n',f_opt_reg,stats_homotopy.comp_res,stats_homotopy.n_biactive,stats_homotopy.cpu_time_total,stats_homotopy.success,stats_homotopy.multiplier_based_stationarity)
 fprintf('Active Set \t %2.2e \t %2.2e \t\t %d \t\t\t %2.2f \t\t\t\t %d\t %s\n',f_opt_active_set,stats_active_set.comp_res,stats_active_set.n_biactive,stats_active_set.cpu_time_total,stats_active_set.success,stats_active_set.multiplier_based_stationarity)
 fprintf('-------------------------------------------------------------------------------\n');
-fprintf(' || w_homotopy - w_active_set || = %2.2e \n',norm(w_opt_homotopy-w_opt_active_set));
+fprintf(' || w_homotopy - w_active_set || = %2.2e \n',norm(w_opt_reg-w_opt_active_set));
+
 
-x_opt_active_set = w_opt_active_set(1:n);
-x_opt_homotopy = w_opt_homotopy(1:n);
 %%
-cardinality_active_set = sum(heaviside(abs(x_opt_active_set)-1e-3))
-cardinality_homotopy  = sum(heaviside(abs(x_opt_homotopy)-1e-3))

src/+mpecopt/Options.m

@@ -38,6 +38,7 @@
         max_recovery_iters (1,1) double {mustBeReal, mustBeInteger} = 12; % If the current BNLP is infeasible, try to solve a tighet relaxation of the MPEC for a better feasible BNLP guess;
         bnlp_projection_strategy (1,1) BNLPProjectionStrategy = BNLPProjectionStrategy.LPEC; % chose how to project x_k(tau) onto complementarity set
         project_guess_to_bounds (1,1) logical = false; % project point to comps and simple bounds (might be still infeasible for general constraints)
+        warm_start_lpec_phase_i (1,1) logical = false; % if true, use latest lpec solution as initial guess (it may be infeasible)
 
 
         relax_and_project_iters (1,1) double {mustBeReal, mustBeInteger} = 1;
@@ -69,6 +70,7 @@
         accept_last_inner_iter(1,1) logical = true; % if max number of inner iterations reached, accept step even if not sufficient decrees      
         smaller_tr_in_phase_ii (1,1) logical = false; % TODO: align with delta tr abooveonce feasiblity is reached, solve lpecs with a very small tr.
         reset_TR_radius (1,1) logical = true;
+        warm_start_lpec_phase_ii (1,1) logical = false; % if true, use latest lpec solution as initial guess (it may be infeasible)
 
         % ----- LPEC solver settings ----
         consider_all_complementarities_in_lpec(1,1) logical = true; % if true, take all comps into the lpec, if false take only those active at x^k

src/+mpecopt/Solver.m

@@ -148,8 +148,10 @@
                             stats.iter.cpu_time_lpec_preparation_iter = [stats.iter.cpu_time_lpec_preparation_iter;toc(t_lpec_preparation_iter)];
                             %  Initial guess and TR for the LPEC
                             y_lpec_k_previous = y_lpec_k_l; % to keep track of active set chnages
-                            lpec.d_lpec = d_lpec_k_l;
-                            lpec.y_lpec = y_lpec_k_l; % inital guess for bin. variablels.
+                            if opts.warm_start_lpec_phase_i
+                                lpec.d_lpec = d_lpec_k_l;
+                                lpec.y_lpec = y_lpec_k_l; % inital guess for bin. variablels.
+                            end
 
                             rho_TR_k_l_lb = opts.realx_and_project_scale_factor_rho_tr*max(abs(min(x_k(dims.ind_x1),x_k(dims.ind_x2))));
                             rho_TR_k_l_ub = opts.realx_and_project_scale_factor_rho_tr*max(max(abs(x_k(dims.ind_x1)),abs(x_k(dims.ind_x2))));

src/mpec_optimizer.m

@@ -127,8 +127,10 @@
                     stats.iter.cpu_time_lpec_preparation_iter = [stats.iter.cpu_time_lpec_preparation_iter;toc(t_lpec_preparation_iter)];
                     %  Initial guess and TR for the LPEC
                     y_lpec_k_previous = y_lpec_k_l; % to keep track of active set chnages
-                    lpec.d_lpec = d_lpec_k_l;
-                    lpec.y_lpec = y_lpec_k_l; % inital guess for bin. variablels.
+                    if settings.warm_start_lpec_phase_i
+                        lpec.d_lpec = d_lpec_k_l;
+                        lpec.y_lpec = y_lpec_k_l; % inital guess for bin. variablels.
+                    end
 
                     rho_TR_k_l_lb = settings.realx_and_project_scale_factor_rho_tr*max(abs(min(x_k(dims.ind_x1),x_k(dims.ind_x2))));
                     rho_TR_k_l_ub = settings.realx_and_project_scale_factor_rho_tr*max(max(abs(x_k(dims.ind_x1)),abs(x_k(dims.ind_x2))));

src/mpecopt_phase_ii.m

@@ -55,8 +55,10 @@
     while ~accept_trail_step && l_k <= settings.max_inner_iter && ~stats.stopping_criterion_fullfiled
         % Here one could do a fast_B_stationarity_check
         y_lpec_k_previous = y_lpec_k_l; % to keep track of active set chnages
-        lpec.d_lpec = d_lpec_k_l; % Initial guess and TR for the LPEC
-        lpec.y_lpec = y_lpec_k_l; % inital guess for bin. variablels.
+        if settings.warm_start_lpec_phase_ii
+            lpec.d_lpec = d_lpec_k_l; % Initial guess and TR for the LPEC
+            lpec.y_lpec = y_lpec_k_l; % inital guess for bin. variablels.
+        end
         lpec.rho_TR = rho_TR_k_l; % update trust region
         stats.iter.rho_TR_iter = [stats.iter.rho_TR_iter, rho_TR_k_l]; % store TR radius
         % Solve LPEC