From 33c656fdf53e7198b8a1025516d7a7905e39d740 Mon Sep 17 00:00:00 2001 From: Anton Edvinovich Pozharskiy Date: Fri, 28 Feb 2025 16:51:57 +0100 Subject: [PATCH 1/4] don't warmstart the lpec with prev --- src/mpecopt_phase_ii.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mpecopt_phase_ii.m b/src/mpecopt_phase_ii.m index 251c294..c0f6488 100644 --- a/src/mpecopt_phase_ii.m +++ b/src/mpecopt_phase_ii.m @@ -55,7 +55,7 @@ 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.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. 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 From a3e0c43ce5d4360b4f0d333f5045a5bc53b17d48 Mon Sep 17 00:00:00 2001 From: "armin.nurkanovic" Date: Thu, 6 Mar 2025 11:36:24 +0100 Subject: [PATCH 2/4] addapt cardinality example --- examples/cardinality_optimization.m | 44 ++++++++++++++++++----------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/examples/cardinality_optimization.m b/examples/cardinality_optimization.m index 4a92034..2edc71c 100644 --- a/examples/cardinality_optimization.m +++ b/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)) From baef82182eca81f8200663a0acb7460d4c394e99 Mon Sep 17 00:00:00 2001 From: "armin.nurkanovic" Date: Thu, 6 Mar 2025 12:54:02 +0100 Subject: [PATCH 3/4] make warm starts of lpec optional --- src/+mpecopt/Options.m | 2 ++ src/+mpecopt/Solver.m | 6 ++++-- src/mpecopt_phase_ii.m | 6 ++++-- 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/+mpecopt/Options.m b/src/+mpecopt/Options.m index b22a900..f2dcc8d 100644 --- a/src/+mpecopt/Options.m +++ b/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 diff --git a/src/+mpecopt/Solver.m b/src/+mpecopt/Solver.m index 5b61ce9..740c48b 100644 --- a/src/+mpecopt/Solver.m +++ b/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)))); diff --git a/src/mpecopt_phase_ii.m b/src/mpecopt_phase_ii.m index c0f6488..cabd9e0 100644 --- a/src/mpecopt_phase_ii.m +++ b/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 From b6569e3fa4b9489dd9f7ebeb9ed9e8800bdc3e00 Mon Sep 17 00:00:00 2001 From: "armin.nurkanovic" Date: Thu, 6 Mar 2025 13:22:42 +0100 Subject: [PATCH 4/4] warm start option in classic mpec_optimizer --- src/mpec_optimizer.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/mpec_optimizer.m b/src/mpec_optimizer.m index 4015ee7..8956036 100644 --- a/src/mpec_optimizer.m +++ b/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))));