SOSTOOLS icon indicating copy to clipboard operation
SOSTOOLS copied to clipboard

Removing initial and final comma to allow for multi-objective optimization

Open rayomaz opened this issue 1 year ago • 4 comments

The initial and final comma in the objvartable of sossetobj.m causes issues. To allow for multi-objective optimization, remove these commas.

rayomaz avatar Aug 06 '24 21:08 rayomaz

Thank you for getting back to this!

Let me re-phrase what I am trying to do here. In the context of the objective, I was referring to multiple optimization variables, that are affine in a single objective.

Simply, I am trying to solve the following program

objfunc = gamsym + betasym;
prog = sossetobj(prog, objfunc);

SOSTOOLS tries to list the optimization variables as a string '[gamsym, betasym]'. Instead, what is returned is ',[gamsym, betasym],'. In other words, a comma is added to the front and end of the string, which throws an error later down the line. From what I can tell, this only happens when having multiple optimization variables in the objective.

Hope that clarifies things.

rayomaz avatar Sep 06 '24 03:09 rayomaz

Thank you for the clarification. However, I don't seem to be able to replicate your issue. I tested some (simple) examples with objective functions in multiple decision variables, but I'm not encountering any error. Also, from what I can tell, although the commas are somewhat awkward, they are explicitly used in the code, and the function doesn't work properly without them. For example, the length of idxcomma1 becomes shorter when you remove the first and last comma, so that the first and last decision variable get skipped in the for-loop starting in line 79.

Of course, you are encountering an error, so perhaps I am missing something. Would it be possible to get just a simple example of a program for which you are encountering an error? If I try just a very rudimentary program

syms x1 x2 gamsym betasym;
prog = sosprogram([x1,x2]);
prog = sosdecvar(prog,[gamsym;betasym]);
prog = sosineq(prog,(x1^2-gamsym));
prog = sosineq(prog,((x2^2+1)-betasym));
prog = sossetobj(prog, -gamsym-betasym);
prog = sossolve(prog);

everything seems to work fine for me without your proposed update.

djagt avatar Sep 06 '24 21:09 djagt

Thanks again.

Here is the program, taken and modified from a paper that synthesizes a stochastic SOS barrier. The code fails at the last line, due to the earlier defined problem with comma based separation programming.

clc; clear; close;

%% Input
alpha_ = 1;   % Range of alpha values for AB(x) =< -alpha*B(x) + Beta
gx = 0.01;    % Range of noise values

N = 1;        % Time horizon (max time)
deg = 2;      % Degree of Polynomials we are searching

%% Function
% Declare symbolic variables
syms z x1 x2 betasym gamsym real;

vars = [x1, x2, betasym, gamsym];
% "z" is the noise variable for which we compute the moments
EXP = 0;

betacoeff = .5;
pJ = .5;
pA = .95;

% Setup Solver
solver_opt.solver = 'mosek';
fx = [  betacoeff*x1; ...
        pJ*x1 + pA*x2; ...
        ];            % Diffusion terms
% Define program
prog = sosprogram(vars);

Zmon = monomials([x1, x2], [0:deg]);
[prog, B] = sospolyvar(prog, Zmon, 'wscoeff');
[prog, sig_u] = sospolyvar(prog, Zmon);
[prog, sig_x] = sospolyvar(prog, Zmon);
[prog, sig_o] = sospolyvar(prog, Zmon);

% Define Inequalities
prog = sosineq(prog, betasym);
prog = sosineq(prog, sig_u);
prog = sosineq(prog, sig_x);
prog = sosineq(prog, sig_o);
prog = sosineq(prog, B);

prog = sosineq(prog, gamsym);
prog = sosineq(prog, 1 - gamsym - 1e-6);
prog = sosineq(prog, B - sig_u*(x1^2 + x2^2 - 2^2 - 1e-3) - 1);
prog = sosineq(prog, -B - sig_o*(1.5^2 - x1^2 - x2^2) + gamsym);

stdvar = gx;
x1 = fx(1);
x2 = fx(2) + z;

Bsub = expand(subs(B));
clear x1 x2;
syms x1 x2 real;
termlist = children(Bsub);

% The following equations are for computer the Expected value
for ii = 1:length(termlist)
    zcount = 0;
    x1count = 0;
    x2count = 0;

    % Initialize expectation as zero value
    EXPz = 0;
    
    factored = cell2sym(termlist(ii));

    factoredterm = factor(factored);
    % Count power of each variable in monomial
    for jj = 1:length(factoredterm)
        
        % Count the power of the "noise" term
        if isequaln(factoredterm(jj),z)
            zcount = zcount + 1;
        end
        
        % Count nubmer of states
        if isequaln(factoredterm(jj),x1)
            x1count = x1count + 1;
        end
        
        if isequaln(factoredterm(jj),x2)
            x2count = x2count + 1;
        end
    end
    
    % Check to see if symbolic term has no noise variables
    % May or may not have state variables
    if zcount == 0
        EXPz = EXPz + factored;
    end
    
    % Check if there are "noise" terms but not x variables
    if zcount > 0 && x1count == 0 && x2count == 0
        if mod(zcount,2) == 1
            EXPz = 0;
        elseif mod(zcount,2) == 0
            EXPz = prod(factoredterm(find(factoredterm~=z)))*prod([1:2:zcount])*stdvar^zcount;
        end
    end
    
    % Check if there are "noise" terms and x variables
    if zcount > 0 && (x1count > 0 || x2count > 0)
        if mod(zcount,2) == 1
            EXPz = 0;
        elseif mod(zcount,2) == 0
            EXPz = prod(factoredterm(find(factoredterm~=z)))*prod([1:2:zcount])*stdvar^zcount;
        end
    end
    
    EXP = EXP + EXPz;
end

%   Expectation evolution
prog = sosineq(prog, -EXP + B + betasym - sig_x*(2^2 - x1^2 - x2^2));

objfunc = gamsym + betasym;
prog = sossetobj(prog, objfunc);

rayomaz avatar Sep 08 '24 02:09 rayomaz

Thanks for sharing that, it's very helpful. From what I can tell, the issue is that the decision variables betasym and gamsym are declared as independent variables in the optimization program. Replacing the lines

vars = [x1, x2, betasym, gamsym];
prog = sosprogram(vars);

with

vars = [x1, x2];
decvars = [betasym, gamsym];
prog = sosprogram(vars);
prog = sosdecvar(prog,decvars);

the code runs fine for me. In fact, you should be able to declare just

vars = [x1, x2];
decvars = [betasym, gamsym];
prog = sosprogram(vars,decvars);

to initialize the program with the desired independent and decision variables. However, it seems that that functionality does not work as desired, so I will be pushing a commit to the SOSTOOLS400 branch to fix that. Thanks for bringing my attention to this!

Does the code run for you if you declare the decision variables explicitly as decision variables? Otherwise, please let me know, and we'll figure out what's going on.

djagt avatar Sep 08 '24 08:09 djagt

I have worked around this issue in a different way. Thank you for looking into this!

rayomaz avatar Nov 15 '24 18:11 rayomaz