# Optimization for a multi-summation problem

Dears

I want to use CVX to solve an optimization problem that consist of nested summations
I wrote a code but the solver says that the solution is infeasible. I want to make sure that I did it right or not. “sorry for not knowing to to make the whole code appear like a code in your website”

My code is

``````	% Optimal Price Regulation Based Algorithm for 500 EVs
%% Updating the SOC and CR at every 5 minute and PD set to zero if the SOC
% is greater than the 0.99Mci and running the simulation for the whole 3
% months from 21 July, 2010 to 20 Oct, 2010

%%
clc
clear

%% Initializing Variables

% NEV  = 500;
NEV  = 10;
Ef   = 0.9;
days = 2;
Beta=0.05;
hours = 4;
weeks=3;
years=12;
BatC=200; %200\$/KWH so the total cost is MC*Batc
discount_rate=0.05;
Meter_cost=36.1; %29 euro
Comm_cost=88.4; %71 euro
Bidirecional_cost=186.7; %150 euro
weighted_factor=[4 4 4]; %3 months, 4 weeks per month
ExD = zeros(hours,days,weeks);
ExU = zeros(hours,days,weeks);
ExR = zeros(hours,days,weeks); %responsive reserve
Pregup=zeros(days,hours,weeks);
Pregdw=zeros(days,hours,weeks);
Prespres=zeros(days,hours,weeks);
POP  = zeros(NEV,hours,days,weeks);
MnAP = zeros(NEV,hours,days,weeks);
MxAP = zeros(NEV,hours,days,weeks);
RsRP = zeros(NEV,hours,days,weeks); %responsive reserve
Deg = zeros(NEV,hours,days,weeks);
% EPD   = zeros(NEV,hours,days);
% ReguppD = zeros(hours,days,weeks);
% RegdwpD = zeros(hours,days);
% ResrespD = zeros(hours,days);
PriceD = zeros(hours,days,weeks);
FPriceD = zeros(hours,days,weeks);
DC=zeros(NEV,1);
ReguppW=zeros(hours*days,weeks);
RegdwpW=zeros(hours*days,weeks);
ResrespW=zeros(hours*days,weeks);

ExU_W=zeros(hours*days,weeks);
ExD_W=zeros(hours*days,weeks);
ExR_W=zeros(hours*days,weeks);

PriceD_W=zeros(hours*days,weeks);
FPriceD_W=zeros(hours*days,weeks);

infeasiable_Days=zeros(days,1);
failed_Days=zeros(days,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SOCI = [1.187409836
1.484262295
0.593704918
1.042622951
0.890557377
1.088459016
0.692655738
1.385311475
2.374819672
2.374819672
];
MP = [3.3
3.3
3.3
3.3
3.3
3.3
3.3
3.3
3.3
3.3
];
Mci = [21.6
21.6
21.6
21.6
21.6
21.6
21.6
21.6
21.6
21.6
];
Av=[1	1	0	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	0	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	1	1	0	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	1	0	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	1	1	0	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	1	0	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	0	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	1	1	0	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1
1	1	0	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1	1
1	1	0	1	1	1	1	1	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1
];
% % Reading the Morning Trip Time of EVs
Mttrip = [3
2
4
3
4
3
2
4
3
3
];
% % Reading the Evening Trip Time of EVs
Ettrip = [11
11
11
11
11
11
11
13
11
12
];
% % Reading the Trip of EVs of Morning trip (Reduction in SOC)
trip_1 = [2.57536
5.47264
0.035555552
1.93152
0.96576
0.64384
1.28768
2.25344
2.57536
3.54112
];
% % Reading the Trip of EVs of Evening trip (Reduction in SOC)
trip_2= [2.66284221
5.56012221
0.123037762
2.01900221
1.05324221
0.73132221
1.37516221
2.34092221
2.66284221
3.62860221
];
% % Reading the compensation factor of the ith EV to account for unplanned
% departures
total_comp = [1	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1	1	1
1	1	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1	1	1
1	1.002506266	1.002506266	1.002506266	1.002506266	1.002506266	1.002506266	1.002506266	1.002506266	1.002506266	1.002506266	1	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1	1	1
1	1	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1	1	1
1	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1.002277904	1	1.006289308	1.006289308	1.006289308	1.006289308	1.006289308	1.006289308	1.006289308	1.006289308	1	1	1
1	1	1	1	1.0041841	1.0041841	1.0041841	1.0041841	1.0041841	1.0041841	1	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1	1	1
1	1	1	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1.005586592	1	1	1
1	1	1	1.003584229	1.003584229	1.003584229	1.003584229	1.003584229	1.003584229	1.003584229	1	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1	1	1
1	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1.002785515	1	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1	1	1
1	1	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1.003134796	1	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1.005025126	1	1	1
];
% % Reading the Expected percentage of the EVs remaining to perform V2G at
% hour t
EvPer = [1	0.989671717	0.980991162	0.974294733	0.970128066	0.970128066	0.970128066	0.970128066	0.970128066	0.970128066	0.986546717	0.972727273	0.952777778	0.946527778	0.946527778	0.946527778	0.946527778	0.946527778	0.946527778	0.946527778	0.946527778	1	1	1
];

%% Price Data from XLS

FPrice = [25.34844871	26.0056763	32.12181122	35.1589504	41.49888019	40.37125422	48.02618913	99.95810515	170.1224283	114.3232142	81.99950949	67.0534449	45.06777494	38.82448189	36.39142323	40.15932157	38.78563793	37.70701738	43.54440679	33.91649244	28.2987497	22.82627733	23.75272459	24.66849407
];
Regupp = [3.67	4.69	6.03	9.85	9.12	7.03	10	13.67	17.03	20.95	20	20	9.67	10	12	5	7.71	5.67	15	4	3	3	3	3
];
Regdwp = [4.37	3.81	4.77	4.62	3.81	5.47	6.68	7.07	7.5	10	10	9.38	6.14	4.65	3.49	4.99	6.3	4.27	8	5.03	4	7.03	6.5	5.03
];

Resresp = [4	3.77	5.47	5.48	5	8	10	9.37	15	24.56	25.64	26.47	10.03	9.49	12	5	7.69	5.67	6.32	3.32	2.69	2.69	2.69	2.76
];

%% ExU and ExD Data from XLS

ExUA = [0.015632083	0.045573417	0.056452813	0.114570522	0.31130082	0.032311417	0.186061917	0.241440833	0.294051667	0.378971333	0.014297208	0.164136661	0.452626392	0.26786129	0.0804151	0.042010858	0.065476775	0.002456933	0.045135388	0.35239	0.853553333	0.217934583	0.130421167	0.13385645
];
ExDA = [0.357682712	0.169056759	0.252310833	0.036405331	0.190553	0.229229617	0.236906592	0.168317187	0.310912917	0.20496615	0.3485209	0.337526083	0.009581686	0.192533192	0.400753607	0.430311917	0.36266675	0.434689667	0.223046766	0.03743375	0	0.526098833	0.065650442	0.038289736
];

ExRA = [0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic

t=hours;
for i=1:NEV
Investment_cost(i)=Meter_cost+Comm_cost+(Bidirecional_cost*MP(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cvx_begin %quiet %quiet for not producing any screen outputs
%in variables we put decision and axuillary variables that are related to
%decision variables
variables POP(NEV,t,days,weeks) MxAP(NEV,t,days,weeks) MnAP(NEV,t,days,weeks) RsRP(NEV,t,days,weeks) Deg(NEV,t,days,weeks)
expressions Inn(t,days,weeks) running_cost(t,days,weeks)  Scaled_week_profit(weeks) Yearly_profit(years) Total_profit

for week=1:weeks
week
%loading data for the whole week (the week here is 2 days, each is 4
%hours for simplicity)
ReguppW(:,week) = Regupp((week-1)*days*hours+1:week*days*hours);
RegdwpW(:,week) = Regdwp((week-1)*days*hours+1:week*days*hours);
ResrespW(:,week) = Resresp((week-1)*days*hours+1:week*days*hours);

ExU_W(:,week) = ExUA((week-1)*days*hours+1:week*days*hours);
ExD_W(:,week) = ExDA((week-1)*days*hours+1:week*days*hours);
ExR_W(:,week) = ExRA((week-1)*days*hours+1:week*days*hours);

FPriceD_W(:,week) = FPrice((week-1)*days*hours+1:week*days*hours);

for day = 1:days

day

Pregup(day,:,week) = ReguppW((day-1)*hours+1 : day*hours,week)'*1e-3;

% RegdwpD(:,day,week) = RegdwpW((day-1)*hours+1 : day*hours,week);

Pregdw(day,:,week) = RegdwpW((day-1)*hours+1 : day*hours,week)'*1e-3;

% ResrespD(:,day,week) = ResrespW((day-1)*hours+1 : day*hours,week);
%Pregup = ReguppD(srt_hrs:end_hrs,day)'*1e-3;   % *1e-3 for KW to MW for 9
%hours
Prespres(day,:,week) = ResrespW((day-1)*hours+1 : day*hours,week)'*1e-3;   % *1e-3 for KW to MW for 24 hours Responsive reserve price

ExU(:,day,week) = ExU_W((day-1)*hours+1 : day*hours,week); %for the whole day
ExD(:,day,week) = ExD_W((day-1)*hours+1 : day*hours,week);
ExR(:,day,week) = ExR_W((day-1)*hours+1 : day*hours,week);

FPriceD(:,day,week) = FPriceD_W((day-1)*hours+1:day*hours,week)*1e-3; %for the whole day

for i=1:NEV
DC(i,1)=0.042*((BatC*Mci(i))/5000)+ ((1-Ef^2)/Ef)*Beta;
end

% the for loop was done to apply the constraints for each time
for k = 1:t %number of hours

Comp(:,k) = (total_comp(:,k));
k

Inn(k,day,week) == ((Pregup(day,k,week)*(sum(MnAP(:,k,day,week))) + Pregdw(day,k,week)*(sum(MxAP(:,k,day,week))) +Prespres(day,k,week)*(sum(RsRP(:,k,day,week))))*EvPer(k))...
+ Beta*(Av(:,k)'*((MxAP(:,k,day,week)*ExD(k,day,week) + POP(:,k,day,week) - MnAP(:,k,day,week)*ExU(k,day,week) - RsRP(:,k,day,week)*ExR(k,day,week))*EvPer(k)));

running_cost(k,day,week)== (Av(:,k)'*((MxAP(:,k,day,week)*ExD(k,day,week) + POP(:,k,day,week) - MnAP(:,k,day,week)*ExU(k,day,week) - RsRP(:,k,day,week)*ExR(k,day,week))*(FPriceD(k,day,week))...
*EvPer(k)))+ (sum(Deg(:,k,day,week)));

end

end

Scaled_week_profit(week)=weighted_factor(week)*(( sum(sum(Inn(:,:,week)))- sum(sum(running_cost(:,:,week))) ));
end

for year=1:years
Yearly_profit(year)=(1/(1+(discount_rate^year)))*((sum(Scaled_week_profit)));
end

Total_profit=sum(Yearly_profit)-sum(Investment_cost);

maximize Total_profit
subject to

for week=1:weeks
week
for day =1:days
day
for k=1:t
k
for i = 1:NEV

if k<Mttrip(i)
Trip(i)=0;
elseif k>=Mttrip(i) && k<Ettrip(i)
Trip(i)=trip_1(i);
elseif k>=Ettrip(i)
Trip(i)=trip_1(i)+trip_2(i);
end

sum(Av(i,1:k).*(MxAP(i,1:k,day,week).*ExD(1:k,day,week)' + POP(i,1:k,day,week) - MnAP(i,1:k,day,week).*ExU(1:k,day,week)' - RsRP(i,1:k,day,week).*ExR(1:k,day,week)').*Comp(i,1:k) + ...
((Deg(i,1:k,day,week)/DC(i,1))*((1-Ef^2)/Ef)))*Ef + SOCI(i) -Trip(i) <= Mci(i);

sum(Av(i,1:k).*(MxAP(i,1:k,day,week).*ExD(1:k,day,week)' + POP(i,1:k,day,week) - MnAP(i,1:k,day,week).*ExU(1:k,day,week)' - RsRP(i,1:k,day,week).*ExR(1:k,day,week)').*Comp(i,1:k) + ...
((Deg(i,1:k,day,week)/DC(i,1))*((1-Ef^2)/Ef)))*Ef + SOCI(i) -Trip(i) >= 0.1*Mci(i); %avoiding deep discharging

if k==t
sum(Av(i,1:k).*(MxAP(i,1:k,day,week).*ExD(1:k,day,week)' + POP(i,1:k,day,week) - MnAP(i,1:k,day,week).*ExU(1:k,day,week)' - RsRP(i,1:k,day,week).*ExR(1:k,day,week)').*Comp(i,1:k) + ...
((Deg(i,1:k,day,week)/DC(i,1))*((1-Ef^2)/Ef)))*Ef + SOCI(i) -Trip(i)>= 0.99*Mci(i);
sum(Av(i,1:k).*(MxAP(i,1:k,day,week).*ExD(1:k,day,week)' + POP(i,1:k,day,week) - MnAP(i,1:k,day,week).*ExU(1:k,day,week)' - RsRP(i,1:k,day,week).*ExR(1:k,day,week)').*Comp(i,1:k) + ...
((Deg(i,1:k,day,week)/DC(i,1))*((1-Ef^2)/Ef)))*Ef + SOCI(i) -Trip(i)<= Mci(i);
end

(MxAP(i,k,day,week) + POP(i,k,day,week))*Comp(i,k) <= (Av(i,k)*MP(i));
MnAP(i,k,day,week)*Av(i,k) <= POP(i,k,day,week) + MP(i)*Av(i,k);
RsRP(i,k,day,week) <= POP(i,k,day,week) + MP(i)*Av(i,k) - MnAP(i,k,day,week);
MxAP(i,k,day,week) >= 0;
MnAP(i,k,day,week) >= 0;
RsRP(i,k,day,week) >= 0;
POP(i,k,day,week) >= -MP(i)*Av(i,k);
Deg(i,k,day,week)>=0;

Deg(i,k,day,week)>= -(DC(i,1)*min(0,Av(i,k)*(POP(i,k,day,week) - MnAP(i,k,day,week)*ExU(k,day,week)- RsRP(i,k,day,week)*ExR(k,day,week)))*(Comp(i,k)/Ef));
end
end
end
end

cvx_end
``````

The output using SDPT3 solver is

## Calling SDPT3 4.0: 3168 variables, 1380 equality constraints For improved efficiency, SDPT3 is solving the dual problem.

num. of constraints = 1380
dim. of linear var = 3120
dim. of free var = 48 *** convert ublk to lblk

SDPT3: Infeasible path-following algorithms

## number of iterations = 13 residual of dual infeasibility certificate X = 1.35e-014 reldist to infeas. <= 1.27e-015 Total CPU time (secs) = 0.30 CPU time per iteration = 0.02 termination code = 2 DIMACS: 1.1e-004 0.0e+000 7.4e+000 0.0e+000 -1.0e+000 1.1e-001

Status: Infeasible
Optimal value (cvx_optval): -Inf

## and using sedumi solver Calling SeDuMi 1.34cvx: 3168 variables, 1380 equality constraints For improved efficiency, SeDuMi is solving the dual problem.

SeDuMi 1.34cvx by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Split 48 free variables
eqs m = 1380, order n = 3217, dim = 3217, blocks = 1
nnz(A) = 13236 + 0, nnz(ADA) = 55528, nnz(L) = 44786
it : by gap delta rate t/tP t/tD* feas cg cg prec
0 : 3.67E+002 0.000
1 : 0.00E+000 2.33E+002 0.000 0.6357 0.9000 0.9000 7.01 1 1 2.7E+000
2 : 0.00E+000 1.33E+002 0.000 0.5708 0.9000 0.9000 2.02 1 1 1.9E+000
3 : 0.00E+000 5.35E+001 0.000 0.4018 0.9000 0.9000 0.44 1 1 2.2E+000
4 : 0.00E+000 1.42E+001 0.000 0.2660 0.9000 0.9000 -0.51 1 1 4.9E+000
5 : 0.00E+000 9.87E-001 0.000 0.0694 0.9900 0.9900 -0.88 1 1 1.0E+000
6 : 0.00E+000 2.77E-004 0.000 0.0003 0.9999 0.9999 -0.99 1 1
Dual infeasible, primal improving direction found.
iter seconds |Ax| [Ay]_+ |x| |y|
6 0.9 1.3e-009 2.7e-013 3.2e+000 1.1e-012

## Detailed timing (sec) Pre IPM Post 1.608E+000 2.221E+000 8.500E-002 Max-norms: ||b||=0, ||c|| = 3.451492e+001, Cholesky |add|=0, |skip| = 0, ||L.L|| = 13.036.

Status: Infeasible
Optimal value (cvx_optval): -Inf

You haven’t provided a reproducible example, complete with input data.

Thanks for your reply…there are excel sheets that contains the inputs…however I will put the whole code for you…I really appreciate your help

As of now, your data (EV_Data_100) doesn’t appear to be accessible.

EV_Data_100 are data from an excel sheet…I want to ask about my writing of the problem formulation. I want to maximize the total profit which consists of nested summations(years, weeks,days, no of EVs), then after “subject to” I wrote the problem constraints…I used also “expressions”

so L don’t know whether this is the way to do things or I have a syntax error???

Manually checking your entire formulation and code requires a bit more energy than I am prepared to spend, and I don’t even know what your model is supposed to be. I can not reproduce your results since I don’t have your input data. Maybe someone else will be able to help you in these circumstances.

Ok. I will prepare a code with an input sample and does not require excel sheets…do you want me to put it in a new question???

Just put it here. I suggest you edit your post in this thread. While you’re at it, try a couple of different solvers (cvx_solver) and post the output from running your model with each solver.

I inserted a new code with ready inputs and with the output of Sedumi and SDPT3 solvers…Thanks for your help

I have reproduced your results. I’ll let mcg weigh in with a better assessment than I can make.

Thank you and I hope to hear from mcg soon…

Dears… I want to ask if there is any update regarding my problem?

Dears… I want to ask if there is any update regarding my problem?

This is not a real answer, but is being posted to draw mcg’s attention to the comment stream. Thanks.

I see no reason to distrust CVX here. Both solvers say the model is infeasible. I can’t tell you what’s wrong with the model, because it’s not my model. You need help not from CVX people, but from people who understand your model.

Got it, but there is really nothing to add…

I mean from the CVX point of view…the way that the nested summations were handled is correct??? is CVX capable of solving 4-dimensional decision variable problem??? I want to make sure about this so I can remove a syntax error or the incapability of CVX in solving this problem??

You have declared 4 dimensional decision variables for which CVX has solved (even if the result is infeasible). I have successfully solved my own 4 dimensional problems in CVX. Your model does not generate any CVX error messages. (continued)

However, as to whether you have correctly implemented your intended model, you hopefully know better than us what you intended, and if you don’t, then you need help from those who do (as mcg was saying above).