Wrong solution using bnb, fmincon can not solve non-linear optimization

[Sunandan Adhikary]

I formulated the following optimization problem. To my understanding, it is a non-linear one ( hope a convex one, if non-convex, I have solver for that too). But YALMIP recognizes it as integer branch & bound. Tried with fmincon, gurobi both, but still I guess since they do not get any feasible solution, YALMIP invokes bnb (even bmibnb does not get invoked if I choose via sdpsettings!).

THE FORMULATION:

sdpsettings('cachesolvers',1'verbose',1,'solver','fmnicon');

binvar a;

binvar b;

sdpvar e;

K1 = sdpvar(1,2);

x = sdpvar(4,1);

A = 0.1*eye(2,2);

B = [0.05;4];

C = [0 1];

K = [0.2  0.09];

L = [-0.04;0.45];

lim0 = [1;2;1;2];

lim1 = 3;

lim2 = 1;

in_reg = 0.1;

AA= [A zeros(size(A));L*C (A-L*C)] ;

BB= [B;B];

C1 = [C zeros(1,size(A,1))];

C2 = [zeros(1,size(A,1)) C];

KK= -[zeros(size(K)) K+K1];

AA_cl= AA+BB*KK;

initlimlow = [a*(-lim0(1)) + (1-a)*(lim0(1)-e); %a=0=> lim0(1)-e<=lim0(1)

                b*(-lim0(2)) + (1-b)*(lim0(2)-e);

                -lim0(1);

                -lim0(2)];            

initlimup = [a*(-lim0(1)+e)+(1-a)*lim0(1);

                b*(-lim0(2)+e)+(1-b)*lim0(2);

                lim0(1);

                lim0(2)];

th = 4.35;

cov = 0.1169;

t=2;

constraints = [initlimlow <= x <= initlimup];

constraints = [constraints, 0 <= e <= 1];

constraints = [constraints, (-1)*sqrt(th*nonatk_cov) <= C1*x - C2*x <= sqrt(th*nonatk_cov)];

constraints = [constraints, (-1)*lim1 <= C1*x <= lim1];

constraints = [constraints, (-1)*lim2 <= KK*x <= lim2];

objective = KK*x 

for i=1:t

    constraints = [constraints, (-1)*sqrt(th*cov) <= C1*power(AA_cl,i)*x - C2*power(AA_cl,i)*x <= sqrt(th*cov)];

    constraints = [constraints, (-1)*lim1 <= C1*power(AA_cl,i)*x <= lim1];

    constraints = [constraints, (-1)*lim2 <= KK*power(AA_cl,i)*x <= lim2];

    constraints = [constraints, (-1)*lim0 <= power(AA_cl,i)*x <= lim0];

    objective = objective +KK*power(AA_cl,i)*x;

end

constraints = [constraints, (-1)*in_reg*lim0 <= power(AA_cl,t)*x <= in_reg*lim0];

% sol = optimize(constraints, objective)

sol = optimize(constraints)

value(K1)

OUTPUT:

Warning: YALMIP has detected that your drive or network is unusually slow.

This causes a severe delay in OPTIMIZE when I try to find available solvers.

To avoid this, use the options CACHESOLVERS in SDPSETTINGS.

See the FAQ for more information.

* Starting YALMIP integer branch & bound.

* Lower solver   : fmincon-standard

* Upper solver   : rounder

* Max time       : 3600

* Max iterations : Inf

 

Warning : The continuous relaxation may be nonconvex. This means 

that the branching process is not guaranteed to find a

globally optimal solution, since the lower bound can be

invalid. Hence, do not trust the bound or the gap...

 Node       Upper       Gap(%)     Lower     Open   Elapsed time

    1 :          Inf      Inf      0.000E+00   1       0.0    Successfully solved  

    2 :          Inf      Inf      0.000E+00   2       0.1    Successfully solved  

    3 :    0.000E+00     0.00      0.000E+00   0       0.1    -> Found improved solution! 

+   3 Finishing.  Cost: 0

sol = 

  struct with fields:

    yalmipversion: '20200930'

    matlabversion: '9.9.0.1592791 (R2020b) Update 5'

       yalmiptime: 2.352955

       solvertime: 0.0815269999999995

             info: 'Successfully solved (BNB)'

          problem: 0

ans =

  Column 1

        0.0634650856801199

  Column 2

        -0.281171683888225

[Johan Löfberg]

First, you call sdpsettings but do not assign it to any variable, nor do you use it in the call to optimize

and it is definitely not convex as you've got loads of general polynomial stuff in the constraints

...and it does not run

Unrecognized function or

variable 'nonatk_cov'.

 

[Sunandan Adhikary]

I am sorry. Please find the proper code below. GUROBI and fmincon does not support polynomial constraints so ran with bmibnb. 

clear all;

clc;

ss= sdpsettings('verbose',1,'solver','bmibnb');

binvar a;

binvar b;

sdpvar e;

K1 = sdpvar(1,2);

x = sdpvar(4,1);

% A = 0.1*eye(2,2);

% B = [0.05;4];

% C = [0 1];

% K = [0.2  0.09];

% L = [-0.04;0.45];

A = [1.0000    0.1000; 0    1.0000];

B = [0.0050; 0.1000];

C = [1 0];

D = [0];

K = [16.0302    5.6622];  % settling time around 10

L = [0.9902; 0.9892];

lim0 = [25;30;25;30];

lim1 = 30;

lim2 = 36;

in_reg = 0.3;

AA= [A zeros(size(A));L*C (A-L*C)] ;

BB= [B;B];

C1 = [C zeros(1,size(A,1))];

C2 = [zeros(1,size(A,1)) C];

KK= -[zeros(size(K)) K+K1];

AA_cl= AA+BB*KK;

initlimlow = [a*(-lim0(1)) + (1-a)*(lim0(1)-e); %a=0=> lim0(1)-e<=lim0(1)

                b*(-lim0(2)) + (1-b)*(lim0(2)-e);

                -lim0(1);

                -lim0(2)];            

initlimup = [a*(-lim0(1)+e)+(1-a)*lim0(1);

                b*(-lim0(2)+e)+(1-b)*lim0(2);

                lim0(1);

                lim0(2)];

th = 4.35;

cov = 0.1169;

t=2;

constraints = [initlimlow <= x <= initlimup];

constraints = [constraints, 0 <= e <= 1];

%% 1st iteration

constraints = [constraints, (-1)*sqrt(th*cov) <= C1*x - C2*x <= sqrt(th*cov)];

constraints = [constraints, (-1)*lim1 <= C1*x <= lim1];

constraints = [constraints, (-1)*lim2 <= KK*x <= lim2];

objective = KK*x 

for i=1:t

    constraints = [constraints, (-1)*sqrt(th*cov) <= C1*power(AA_cl,i)*x - C2*power(AA_cl,i)*x <= sqrt(th*cov)];

    constraints = [constraints, (-1)*lim1 <= C1*power(AA_cl,i)*x <= lim1];

    constraints = [constraints, (-1)*lim2 <= KK*power(AA_cl,i)*x <= lim2];

    constraints = [constraints, (-1)*lim0 <= power(AA_cl,i)*x <= lim0];

    objective = objective +KK*power(AA_cl,i)*x;

end

% performance region criteria

constraints = [constraints, (-1)*in_reg*lim0 <= power(AA_cl,t)*x <= in_reg*lim0];

% sol = optimize(constraints, objective)

sol = optimize(constraints,[],ss)

% vpa(value(K1))

value(K1)

Now getting infeasible solution:

Warning: YALMIP has detected that your drive or network is unusually slow.

This causes a severe delay in OPTIMIZE when I try to find available solvers.

To avoid this, use the options CACHESOLVERS in SDPSETTINGS.

See the FAQ for more information.

* Starting YALMIP global branch & bound.

* Upper solver     : fmincon

* Lower solver     : GUROBI

* LP solver        : GUROBI

* -Extracting bounds from model

* -Perfoming root-node bound propagation

* -Calling upper solver Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  7.890194e-19. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  1.091001e-20. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  5.444503e-20. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  6.679043e-22. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  5.904115e-23. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  5.248643e-21. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  5.993699e-21. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  3.013660e-21. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  5.433079e-22. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  2.647184e-23. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

Warning: Matrix is close to singular

or badly scaled. Results may be

inaccurate. RCOND =  2.647184e-23. 

> In backsolveSys

In solveKKTsystem

In computeTrialStep

In barrier

In fmincon (line 848)

In callfmincon (line 95)

In global_solve_upper (line 84)

In solve_upper_in_node (line 9)

In bmibnb (line 214)

In solvesdp (line 371)

In optimize (line 31)

In test (line 68)

 

(no solution found)

* -Branch-variables : 8

* -More root-node bound-propagation

* -Performing LP-based bound-propagation 

* -And some more root-node bound-propagation

* Timing: 93% spent in upper solver (1 problems solved)

*         0% spent in lower solver (0 problems solved)

*         1% spent in LP-based domain reduction (1 problems solved)

*         0% spent in upper heuristics (0 candidates tried)

sol = 

  struct with fields:

    yalmipversion: '20200930'

    matlabversion: '9.9.0.1592791 (R2020b) Update 5'

       yalmiptime: 2.281057

       solvertime: 6.061407

             info: 'Infeasible problem (BMIBNB)'

          problem: 1

ans =

   NaN   NaN

Should I use some other solver or there is something wrong with the problem formulation.

[Johan Löfberg]

Already this is infeasible

constraints = [initlimlow <= x <= initlimup];

constraints = [constraints, 0 <= e <= 1];

constraints = [constraints, (-1)*in_reg*lim0 <= power(AA_cl,t)*x <= in_reg*lim0];

sol = optimize(constraints,[],ss)

Pretty far from feasible also, as a slack shows that the last one has to be massively violated to be feasible

sdpvar q

constraints = [initlimlow <= x <= initlimup];

constraints = [constraints, 0 <= e <= 1];

constraints = [constraints, (-1)*in_reg*lim0-q <= power(AA_cl,t)*x <= in_reg*lim0+q];

sol = optimize(constraints,q,ss)

>> value(q)

ans =

   16.2120

[Sunandan Adhikary]

You mean to say since there is a significant violation that can be seen from the value of q, the problem is infeasible practically. But I was increasing the number of iterations (i.e. the value of t) in such case. But it still does not satisfy all the constraints even if it provides some solution. So I was wondering if I am using the proper solver and modeling the problem correctly.

[Johan Löfberg]

It cannot become feasible by adding more constraints, and don't understand what you mean with "does not satisfy all the constraints even if it provides some solution"

[Ipsita Koley]

He tries to mean by "does not satisfy all the constraints even if it provides some solution" that:

for example, we try to get a value of y, that satisfies the constraint a<x<b where x=f(y)

After running the solver, it gives a value of y such that the simulation shows x=b+e where e>0

In the context of the problem posted by Sunandan, the last constraint (i.e. performance region constraint) does not satisfy.

P.S.: We are working on the same problem.

[Sunandan Adhikary]

Attached simulation result using optimization solution for the similar code but for different system to properly demonstrate our query. Note that the performance criteria (commented) constraint does not get satisfied even after getting successful solution

[Johan Löfberg]

this must be a different problem because it is easily solved, and the computed solution is very feasible i.e. satisfies all constraints with al large margin

>> check(constraints)

 

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|    ID|                            Constraint|   Primal residual|   Dual residual|

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|    #1|     Elementwise inequality (bilinear)|           0.10927|             NaN|

|    #2|                Elementwise inequality|          0.080774|             NaN|

|    #3|                Elementwise inequality|           0.24087|             NaN|

|    #4|                Elementwise inequality|            1.3039|             NaN|

|    #5|     Elementwise inequality (bilinear)|           0.71001|             NaN|

|    #6|                Elementwise inequality|           0.52375|             NaN|

|    #7|     Elementwise inequality (bilinear)|             2.195|             NaN|

|    #8|   Elementwise inequality (polynomial)|            0.8099|             NaN|

|    #9|     Elementwise inequality (bilinear)|           0.96326|             NaN|

|   #10|     Elementwise inequality (bilinear)|           0.22696|             NaN|

|   #11|   Elementwise inequality (polynomial)|            1.5763|             NaN|

|   #12|   Elementwise inequality (polynomial)|           0.75816|             NaN|

|   #13|   Elementwise inequality (polynomial)|            0.9598|             NaN|

|   #14|   Elementwise inequality (polynomial)|           0.63492|             NaN|

|   #15|   Elementwise inequality (polynomial)|            2.4366|             NaN|

|   #16|   Elementwise inequality (polynomial)|           0.81226|             NaN|

|   #17|   Elementwise inequality (polynomial)|           0.98337|             NaN|

|   #18|   Elementwise inequality (polynomial)|          0.083368|             NaN|

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

[Sunandan Adhikary]

Yes, it does get solved but if you notice the data tips in the simulation plot in the attached plot image and the last constraint, they contradict! What could be the reason behind this? The solver we chose does not even run for several iterations either, i.e. seems to easily reach a solution and stop. But when we actually verify that, the last constraint is not obeyed! 

[Johan Löfberg]

I don't know what a data tips in the simulation plot is.

To see if the solver has returned a feasible solution, simply run check(constraints). If those numbers are all non-negative (or at least larger -10^-6 or so as solver have numerical tolerances) then a feasible solution has been computed and delivered. If that does not match with your plots, you have either solved another problem than you thought (i.e. missing or wrong constraints) or you are plotting the wrong thing

[Sunandan Adhikary]

Thank you, professor. I will verify if I am putting the constraint properly. We are new to modeling such complex optimization problems.  Since YALMIP seemed really versatile and so well documented, and our problem seemed so complex to frame and classify, we chose this. Facing infeasibility several times, we wanted to be sure about the modeling and solver choices. 

I tried to remodel the problem by adding more sdp variables. But now adding a few new constraints I face  "total dimension of C should be > length(b)" error. I can not find out sure how is that violated by adding those new constraints! (the model runs on removing those constraints!).

Sorry for the repeated questions and thanks a lot for bearing with me and responding promptly.

[Sunandan Adhikary]

Also, there are warnings from YALMIP optimize function saying "Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.". Is there any generic reason behind this? I mean, should I modify my model somehow in order to get no such warnings and more accurate solutions?

[Johan Löfberg]

That's an error message you would see from SDPT3 when solving too trivial problems

Your problem is not trivial though as you've now gone from a very hard nonconvex nonlinear program  to an extremely absurdly hard nonconvex nonlinear semidefinite program (for which the only available solver is bmibnb, and it typically works when there are a few nicely bounded variables and low polynomial degreem, with SDP supports being extremely experimental, while you have 36 variables to branch in, most of them unbounded). For some reason, the linear relaxations during the branching which SDPT is used for turns out to be badly behaved

In other words, you're not going to solve this model with off-the-shelf software

[Johan Löfberg]

That's not yalmip it's fmincon

It's fmincon saying it has issues solving this problem. Just generic bad numerics which could come from almost anything (bad model, bad geometry of the feasible space, infeasible, unbounded)

[Sunandan Adhikary]

Thank you, professor. Following your deductions and suggestions, I understand that our modeling is somewhat introducing some nasty nonconvexities and nonlinearities while also having unbounded variables. But our objective or the problem at hand is pretty trivial and straightforward. We are just trying to derive robust controller gains to bring a closed-loop system back to its preferable region ASAP from any point in state space staying safe, while the control actions, system outputs are constrained. But the resultant problem turned out to be not so simple! Any suggestions regarding this or similar models that we could follow among YALMIP examples? We will check for unbounded variables and if anything can be simplified in the model following your suggestions. Your comments were really helpful.

[Johan Löfberg]

Trivial to state is not the same thing as trivial or even possible to solve. Find output feedback K such that A+BKC is stable and -1<=K<=1 is trivial to state but an open problem and most likely an intractable optimization problem...

[Sunandan Adhikary]

I understand. Well, thanks again. I will try to check if the model can be improved or simplified.

[Sunandan Adhikary]

I just had one general question regarding such optimization problems. When I define an initial x0 sdpvar and put a constraint with its upper and lower bounds then the optimal or near-optimal solution that the solver can provide satisfies for all x0 in that range right? (obviously considering the residue corresponding to that constraint is positive). Asking this because sometimes I get prompts like "[Your initial point x0 is not between bounds lb and ub; FMINCON shifted x0 to strictly satisfy the bounds". Is it because there exists no feasible solution for any x0 in that bound and any solver is generally allowed to make such choices in such situations?

[Johan Löfberg]

Hard to say, but one thing that can happen is that you are only initializing the variables you have defined, and have no control over variables that YALMIP might be defining internally for modelling purposes, and sometimes those can not (or simply are not) initialized

[Sunandan Adhikary]

Okay. But if YALMIP is using some temporary or internal variables to model the unknown variables before sending it to the solver, the solver must be dealing with those uninitialized variables in a way so that the provided bound constraints against the user defined unknown variables are respected right?

[Johan Löfberg]

As I said, hard to say anything generic as a lot of stuff happens on the way between you and the internals of a solver, and it depends on the model and the solver. As an example (nothing special about it) if you do geomean(x) and assign an initial value to x, it will basically help nothing, as all the lifting variables that are introduced to generate an SOCP model of geomean is a massive mess and no attempt is made to initialize those, and for the solver it is just as hard so it will probably just set all of them to 0 and then do a standard start from infeasible initial, and your initial guess on x will not be feasible in the context of the other varibles

[Sunandan Adhikary]

Okay. So a good way to model my constraints in this case (your SOCP example) would be to directly introduce a bounding constraint on the geomean as well? 
And I also wanted to ask, after the problem is solved, sdpvars are populated with certain values. If those are within the range I provided in constraint, then can it be considered that the given solution is optimal for the range I provided?

[Johan Löfberg]

No, you should just accept that you have no control of stuff that might be happening under the hood. For instance

x = sdpvar(7,1)

optimize(geomean(x)>=1)

will generate a model with 35 variables involved in 7 socp constraints. There is no way for you to do anything outside to guide the initialization of that model (and of course, here it is a problem which is solved using primal-dual solver, so you cannot initialize them at all as you cannot initialize the dual, and no primal-dual solver can efficiently arm-start anyway...)

And regarding you second question, no if you are using a local solver you have absolutely no guarantee that the solver finds a global solution in the interval of interest, or even that it finds a solution)

[Sunandan Adhikary]

Okay, understood. I have noticed YALMIP shows the prompt whenever the mentioned upper bound solver (or lower bound solver) can find a global solution and if not it is handed over to the local solvers I am guessing and can not commit to any global solution?

[Johan Löfberg]

Don't understand what you are saying. If you are using bmibnb it will repeatedly call local upper bound solvers and lower bound solvers (and bound propagation LP solvers), that's the whole idea in order to hone in on the global solution

[Sunandan Adhikary]

Yes, I understand what you are saying. YALMIP prompts these actions as they happen. But I am trying to understand whether I am prompted or somehow can understand that the solution I am getting is a global one or not.

[Johan Löfberg]

You judge it from the gap/termination code. If it terminates with code 0 it means it has terminated with optimality within specified tolerance, if not it will tell you if it terminated due to max iteration/time etc, and you will have extract the gap using savesolveroutput, or simply look at the log to see what the final gap was, and then judge that too decide what you think about the solution

[Sunandan Adhikary]

Understood. Thank you, Professor.

[Sunandan Adhikary]

Professor, taking your suggestion from sdpt3 forum, I reframed my problem to see if I can find out an initial set in state-space that would guarantee a safe trajectory. But mosek/gurobi produces positive primal residues but the solution violates safety constraints! (notable from plot). Mosek frames it as a linear problem and says solved successfully. If I am using bmibnb then getting a better solution then others, but the gap seems to be 0 from the start! So no way of understanding if the solutions are global or which solver is right to use and which one is wrong!!

[Johan Löfberg]

No idea what the optimization actually is supposed to do. You've computed *one* initial state for which all (finite-length) trajectories and control inputs for a linear system dynamics stay in some region

When I run it, it finds that the initial state x0 works

value(x0)

ans =

     0

     0

     0

     0

so what?

The feasible set is hard to visualize as it is 4d, but here is the projection of all first  2 states that are safe

plot(constraints,x0(1:2))

and first 3 can also be projected visually

plot(constraints,x0(1:3))

the simulation you run is not starting from the feasible set, easily seen by the fact that this is infeasible

optimize(constraints + [x0 == x_init])

[Johan Löfberg]

should say "You've computed *one* initial state for which *the* (finite-length) trajectory and control input for a linear system dynamics stay in some region"

[Sunandan Adhikary]

So this formulation is wrong if I want out find out an initial set in state-space that would guarantee a safe trajectory. Any suggestions on the modeling then? And the solver choice?

Using the safer variable I am trying to find out the fraction of the safe intervals the initial state sdpvar x0 should be bounded within so that the trajectories can stay within a safe interval. 

[Johan Löfberg]

The largest possible feasible initial set is simply the polytope you have defined in the contraints. 

x0 = sdpvar(4,1)

constraints = [];

xi=x0;

ui=KK*xi;

for i=1:t

    constraints = [constraints, (-1)*sensor_limit <= C1*xi <= sensor_limit];

    constraints = [constraints, (-1)*sqrt(threshold*nonatk_cov) <= C1*xi - C2*xi <= sqrt(threshold*nonatk_cov)];

    constraints = [constraints, (-1)*actuator_limit <= ui <= actuator_limit];

    constraints = [constraints, lower_safety <= xi <= upper_safety];

    xi= AA_cl*xi;

    ui=KK*xi;

end

% Feasible set is A*x0 <=b

[A,b] = double(polytope(constraints))

What you seem to have modelled previously is the search of the largest possible initial box (ie. an inner approximation of the actual feasible set). You are missing that you are not looking for *a* x0, you are looking for a box such that all trajectories starting in that box are feasible, i.e. a forall operator on x0, i.e. effectively a robust optimization problem. Since the size of the uncertainty set is a decision variable to be maimized, we can formulate the largest possible box could be done with the following robust optimization problem https://yalmip.github.io/parameterizeduncertainty

z = sdpvar(4,1);

center = sdpvar(4,1);

sdpvar scale

x0 = scale*z + center;

constraints = [];

xi=x0;

ui=KK*xi;

for i=1:t

    constraints = [constraints, (-1)*sensor_limit <= C1*xi <= sensor_limit];

    constraints = [constraints, (-1)*sqrt(threshold*nonatk_cov) <= C1*xi - C2*xi <= sqrt(threshold*nonatk_cov)];

    constraints = [constraints, (-1)*actuator_limit <= ui <= actuator_limit];

    constraints = [constraints, lower_safety <= xi <= upper_safety];

    xi= AA_cl*xi;

    ui=KK*xi;

end

optimize([constraints, uncertain(z), -1 <= z <=1],-scale)

value(scale)

value(center)

However this are all re-inventing the wheel, in a bad way as you only look at a finite horizon. You want to compute the largest invariant set, something which is well established, and there is readily available software to do this (MPT for instance)

[Sunandan Adhikary]

Also tried modeling like this (used initial bounds as variables). Would it solve the issue?

[Johan Löfberg]

polytope requires MPT to be installed, but since no projection is involved, it is really just

base = getbase(sdpvar(constraints));

b = base(:,1);

A = -base(:,2:end);

although this will not be in reduced form so there can be redundant facets in the representation

[Sunandan Adhikary]

Let me look into your suggestion. Thank you.

[Sunandan Adhikary]

I will try to frame the problem using polytopic representation from the example link you provided. So, are you suggesting to get this maximized initial box using interval sets (as I am doing now) is not a good way? Polytope is the ideal way to do so?

[Johan Löfberg]

the largest possible axis-aligned feasible cube is much smaller than the set feasible

>> volume(polytope(constraints))

ans =

   9.7182e+05

>> (2*value(scale))^4

ans =

   1.9575e+03

[Sunandan Adhikary]

Understood your modeling. Trying to frame my other safe, robust feedback gain synthesis problem for all the feasible initial states derived from this problem. Thank you very much, professor.

[Johan Löfberg]

finding a state-feedback to find largest possible control invariant set is also a standard problem, although I believe it cannot be done using polytopic theory as the analysis problem, but one has to work with inner approximation ellipsoidal sets instead, which in the end leads to an LMI

[Sunandan Adhikary]

My actual problem (that I started posting about) is actually finding robust feedback gains that would work on top of normal lqr gain to bring a disturbing system starting from for all safe initial states (that is derived modeling the above problem) back to a reference/goal/performance region (which I derive as an invariant set). And I also framed it as a set of similar constraints and LMIs included for switching stability using Lyapunov. Thanks a lot again for all the help. Noted your suggestion.