New cones: DD, SDD, DSOS, SDSOS. #76

renatoloureiro · 2024-07-19T15:04:14Z

This pull request introduces and implements NEW cone options:

Diagonally Dominant (DD), dd
Scaled Diagonally Dominant (SDD), sdd
Diagonally Scaled Sum-of-Squares (DSOS), dsos
Scaled Diagonally Scaled Sum-of-Squares (SDSOS). sdsos

Pseudo example of cones definition while using casos.sdpsol:

opts.Kx = struct('lin', 4, 'psd', 2, 'dd', 3, 'sdd', 5);
opts.Kc = struct('psd', 3, 'dd', 4, 'sdd', 2);

In the case of relaxing a SOS program using casos.sossol:

opts.Kx = struct('lin', 4,  'sdsos', 1);
opts.Kc = struct('dsos', 3);

Note: This pull request answers issue #74

Key contributions include:

DD/SDD Cone Implementation: Added support for DD/SDD cones in decision variables and constraints, along with multiple DD/SDD cones handling and bug fixes.
Additional Cone Options: Added DSOS and SDSOS cones to the polynomial cone class, including methods for their dimensions and mappings.
Enhanced Demos: Updated demos to incorporate DD, SDD, DSOS, and SDSOS cones, showcasing their optimal solutions and integration into the problem-solving process.
(see the attached file debug_dd_sdd.zip that contains additional programs to verify the correct performance of the new cones)
Bug Fixes and Enhancements: Addressed bugs related to the new cones, optimized dimensions, and refined plot outputs.
Documentation and Readme: Updated README to document the new cone options and their functionalities.

Example of new cones:

Polynomial optimization with dsos or sdsos:

x = casos.Indeterminates('x');                % indeterminate variable
f = x^4 + 1*x;                                % some polynomial
g = casos.PS.sym('g');                        % scalar decision variable
sos = struct('x',g,'f',g,'g',f+g);            % define SOS problem
opts = struct('Kc',struct('dsos',1));         % constraint is scalar SOS cone
S = casos.sossol('S','mosek',sos,opts);       % solve by relaxation to SDP
sol = S();                                    % evaluate

Relaxed SDP problem with DD (dd) or SDD (sdd):

ndim = 5;                      % dimension of SDP
x = casadi.SX.sym('x', 4, 1);  % decision variables
A = randn(ndim, ndim);         % create random SDP data
B = randn(ndim, ndim);
A = A+A';
B = B+B';

% define SDP problem
sdp.f = x(1);
sdp.g = eye(ndim) + x(1)*A + x(2)*B;
sdp.g = [sdp.g(:)];
sdp.x = x;

% define cones
opts = struct;
opts.Kx = struct('lin', 4);
opts.Kc = struct('psd', ndim);

S = casos.sdpsol('S','mosek',sdp,opts); % initialize solver
sol = S('lbx',-inf,'ubx',+inf); % solve

1. Current implementation only allows a single DD cone in the free variables, not the constraints 2. Due to the creation of new variables and a change of cones, two additional properties in options class were added: they inherit the upper and lower bounds of the new linear cones that are not requested for the user to define

1. Fix bug related to DD cone of size 1. 2. Ensure that only SX is leveraged

The cones description added were DD, SDD and for the polynomial cones DSOS and SDSOS

Add to sdpopt the cones used at the sdp level, namely, the 'lin', 'psd' and 'dd'

Due to the additional decision variables created from building the DD and constraints, a mapping is defined and is used aposteriori to extract the solution of the original problem formulation from the solved sdp

Make both problems use the same solver for a coherent comparison

1. add comments; 2. add SDD feasible region; 3. add SDD program

This file in this commit is similar to the one related to the DD. currently this file performs SDD to linear as if the constraint was DD.

The function sdd_reduced is added in the file to call the re-building to LP/SOCP.

Check instead if the variable args which is only nonempty if either DD/SDD was present in the original problem formulation.

The solutions obtained are from solving the SDP by either using PSD, DD or SDD cones.

This versio works for a single SDD constraint on the SDP level.

renatoloureiro · 2024-08-07T08:10:57Z

I have fixed the bugs that existed before and created two new demos:

/pendulum: generated the dynamics of a controlled n-link inverted pendulum and its polynomial approximation. A region of attraction program using the SOS, DSOS, or SDSOS is given.
/stats_poly: runs a simple polynomial optimization task for increasing free variable dimension for the different polynomial cones, i.e. SOS, DSOS, SDSOS. Plots the elapsed time and result achieved.

tcunis · 2024-08-14T07:28:03Z

+casos/+package/+solvers/@SdpsolInternal/SdpsolInternal.m

@@ -61,19 +119,18 @@
        % constraint function (vectorized)
        sdp_g = sdp.g(:);

-        if isfield(sdp,'derivatives')
+        if isfield(sdp,'derivatives') && isempty(fieldnames(args))


I don't think the second clause works as intended, for even if there are no DD/SDD constraints, the struct args has the fields dd_lbx etc.

see commit (f245d36) for the requested changes

tcunis

Overall it looks, but there are a few critical issues with the interface of the SDP solver and its embedding into the SOS solver that must be addressed.

+casos/+package/+solvers/@SdpsolInternal/buildproblem.m

+casos/+package/+solvers/@SossdpRelaxation/buildproblem.m

tcunis · 2024-08-14T07:37:15Z

+casos/+package/+solvers/@SossdpRelaxation/buildproblem.m

+if Nsds~=0 || Msds~=0
+sdpopt.Kx.sdd = [Ksdp_x_s(Nds+Ns+1:end); Ksdp_g_s(Mds+Ms+1:end)];  
+end
+


don't you also have to reorder the SDP variables and constraints to ensure PSD < DD < SDD?

+casos/+package/+solvers/@SossdpRelaxation/buildproblem.m

+casos/+package/+solvers/@SdpsolInternal/SdpsolInternal.m

+casos/+package/+solvers/@SdpsolInternal/buildproblem.m

…ous computed derivatives. The flag is one if there was some reduction/relaxation to DD/SDD.

instead of verifying the existence of the field, it verifies if the field is empty.

…deal.

… the SOS level. 1. The output of the SDP is equal to the original setting and not from the generated SDP after DD/SDD reduction/relaxation.

tcunis · 2024-08-21T14:41:00Z

+casos/+package/+solvers/@SdpsolInternal/SdpsolInternal.m

@@ -61,19 +118,18 @@
        % constraint function (vectorized)
        sdp_g = sdp.g(:);

-        if isfield(sdp,'derivatives')
+        if isfield(sdp,'derivatives') && ~flag_reduced 


Could we use the map above to update pre-computed derivatives?

Yes, we can:

The inner if-statement is used cause at the moment I don't find an elegant and short way of doing it.

what about A? Why can't we do the same here?

+casos/+package/+solvers/@SdpsolInternal/buildproblem.m

tcunis · 2024-08-21T14:47:40Z

+casos/+package/+solvers/@SdpsolInternal/buildproblem.m

 % output function
 obj.ghan = casadi.Function('g', ...
            {p x cost lam_a lam_x}, ...
-            {x fx gx lam_x lam_a nan(size(p))}, ...
+            {obj.map.x*x fx obj.map.g*gx obj.map.x*lam_x obj.map.g*lam_a nan(size(p))}, ...


Please check whether the dual variables can simply be mapped by the same mapping as the primal variables/constraints.

tcunis · 2024-08-21T14:57:24Z

+casos/+package/+solvers/@SdpsolInternal/sdd_reduce.m

+opts.Kc = setfield(opts.Kc,'psd', [Mpsd, repmat(2, 1, length(vertcat(M_g{:}))/3+length(vertcat(M_x{:}))/3)]);
+
+% map from new sdp.x to old
+map.x = jacobian(x_original, sdp.x);


given the structure of sdp.x, isn't that map simply blkdiag(speye(Nlin), sparse(0,m), speye(n-Nlin), where m is the number of elements in vertcat(M_g{:}, M_x{:}) and n is the number of elements in the original sdp.x.

tcunis · 2024-08-21T14:59:13Z

+casos/+package/+solvers/@SdpsolInternal/dd_reduce.m

+opts.Kc = setfield(opts.Kc,'lin',Mlin);
+
+% map from new sdp.x to old
+map.x = jacobian(x_original, sdp.x);


given the structure of sdp.x, isn't that map simply blkdiag(speye(Nlin), sparse(0,m), speye(n-Nlin), where m is the number of elements in vertcat(M_g{:}, M_x{:}) and n is the number of elements in the original sdp.x.

tcunis · 2024-08-21T15:00:24Z

It appears that map.x and map.g are the same in dd_reduce and sdd_reduce. Is that correct?

renatoloureiro · 2024-08-21T15:07:13Z

The current implementation uses a mapping from what was previously done. That is, the sdd_reduce receives the SDP already altered in the dd_reduced step. This is done mainly due to the reordering, otherwise I would have a harder time afterwards.

tcunis · 2024-08-22T10:10:44Z

The current implementation uses a mapping from what was previously done. That is, the sdd_reduce receives the SDP already altered in the dd_reduced step. This is done mainly due to the reordering, otherwise I would have a harder time afterwards.

I was referring to the construction of map_DD_2_ORIG_x and map_SDD_2_ORIG_x (or their respective variants for the constraints) individually, viz. (for DD):

sdp.x = [sdp.x(1:Nlin); vertcat(M_g{:}); vertcat(M_x{:}); sdp.x(Nlin+1:end)];

% ...

map.x = jacobian(x_original, sdp.x);
map.g = [speye(length(g_original)), sparse(length(g_original), length(sdp.g)-length(g_original))];

vs. (for SDD):

sdp.x = [sdp.x(1:Nlin); vertcat(M_g{:}); vertcat(M_x{:}); sdp.x(Nlin+1:end)];

% ...

map.x = jacobian(x_original, sdp.x);
map.g = [speye(length(g_original)), sparse(length(g_original), length(sdp.g)-length(g_original))];

It appears that the maps in both cases are equal, which seems counterintuitive.

…approximation)

renatoloureiro and others added 30 commits April 4, 2024 10:20

Enable multiple DD cones in decision variables

0b175f3

DD cone implemented in the constraints

edb52a9

Fix bugs in SdpsolInternal

1674cdd

1. Fix bug related to DD cone of size 1. 2. Ensure that only SX is leveraged

Demo to test DD cone implementation

5728584

Move DD reduction to separate file.

9adcdb9

Support of Casadi v3.6.5: allow for duplicate I/O names in functions.

8a61793

Remove additional options.

9a5c9d9

Reorganized reduction of DD problem.

185c4f8

Set random number generator for DD cone demo.

6972fce

Further improvements of DD reduction code.

f283d7a

Print supported cones.

13c8a0c

Update demo: change output layout

df38a77

Merge: origin/main #73 to solvers DD/SDD development

f6a0080

Add more cones to class Cones

854863b

The cones description added were DD, SDD and for the polynomial cones DSOS and SDSOS

Add DSOS and SDSOS as options to polynomial cones

026e8ce

Get dimensions of dsos and sdsos cones

e93175f

Add to sdpopt the cones used at the sdp level, namely, the 'lin', 'psd' and 'dd'

Enhance demo polyopt with DSOS

2fc3211

Add mapping from solved sdp and the original sdp.

ad8dc92

Due to the additional decision variables created from building the DD and constraints, a mapping is defined and is used aposteriori to extract the solution of the original problem formulation from the solved sdp

Bug fix: same solver for both problem.

7e20511

Make both problems use the same solver for a coherent comparison

Add property 'map' and 'sdd_reduce()' method declaration

291e044

Enhance sdprelaxation demo.

9fb3b74

1. add comments; 2. add SDD feasible region; 3. add SDD program

Base file for SDD to linear/SOCP.

8ea3849

This file in this commit is similar to the one related to the DD. currently this file performs SDD to linear as if the constraint was DD.

Add SDD to matrix_cones property.

eae4c54

The function sdd_reduced is added in the file to call the re-building to LP/SOCP.

Add SDD to demo/polyopt.

0a550dc

Fix string output.

f657357

Delete temporary variable 'flag'.

2e056a1

Check instead if the variable args which is only nonempty if either DD/SDD was present in the original problem formulation.

Enhance demo: add optimal solution obtained to figure.

64bddde

The solutions obtained are from solving the SDP by either using PSD, DD or SDD cones.

SDD reduction: first working iteration.

50260c1

This versio works for a single SDD constraint on the SDP level.

Change solver to Mosek and add figure with graph of resulting functions.

fff126f

Bug fix: correct concatenation.

3410cb4

renatoloureiro marked this pull request as ready for review August 7, 2024 07:55

Small change: pcontour domain and set bisection.

41f3698

tcunis removed the bug Something isn't working label Aug 14, 2024

tcunis reviewed Aug 14, 2024

View reviewed changes

tcunis requested changes Aug 14, 2024

View reviewed changes

renatoloureiro added 8 commits August 14, 2024 13:29

Bug fix: use a boolean flag to decide whether or not to use the previ…

f245d36

…ous computed derivatives. The flag is one if there was some reduction/relaxation to DD/SDD.

Bug fix: incorrect if-statement.

f084026

instead of verifying the existence of the field, it verifies if the field is empty.

Merge branch 'main' of github.com:iFR-ACSO/casos into solvers-DD

b149eaa

Split the variables Ksdp_x_s and Ksdp_g_s using mat2cell, celfun and …

5ae747b

…deal.

Bug fix: mapping from cone reduction only on the SDP level and not on…

3a75017

… the SOS level. 1. The output of the SDP is equal to the original setting and not from the generated SDP after DD/SDD reduction/relaxation.

Enhance: make property as private instead of public.

0a3d60f

Small change: change mdiag to blkdiag.

bc550c8

Update generation of pendulum dynamics to ease computation time.

d823313

tcunis reviewed Aug 21, 2024

View reviewed changes

+casos/+package/+solvers/@SdpsolInternal/buildproblem.m Show resolved Hide resolved

tcunis reviewed Aug 21, 2024

View reviewed changes

renatoloureiro and others added 5 commits August 28, 2024 09:26

Speed up pendulum dynamics generation by using a simplified inverse (…

6af367f

…approximation)

Comment out pcontour plot of sliced solution.

6ccf317

Merge branch 'main' into solvers-DD

0ffd7fb

Enhance: use always sdp.derivates.Jg

6f1bb25

Bug-fix: remove empty fields.

4d9c062

tcunis removed request for Fabian-Geyer and JOlucak January 30, 2025 14:09

New cones: DD, SDD, DSOS, SDSOS. #76

Are you sure you want to change the base?

New cones: DD, SDD, DSOS, SDSOS. #76

Uh oh!

Conversation

renatoloureiro commented Jul 19, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

renatoloureiro commented Aug 7, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tcunis left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tcunis commented Aug 21, 2024

Uh oh!

renatoloureiro commented Aug 21, 2024

Uh oh!

tcunis commented Aug 22, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

renatoloureiro commented Jul 19, 2024 •

edited

Loading

tcunis commented Aug 22, 2024 •

edited

Loading