%% Figure 3  calculation 1 (ST1&ST2)

%%
clear all
close all

load DataNomal


%% setting for calc region/year
set_i_reg=1:10;
set_i_year=1:3;
n_reg=10;
n_year=3;

set_PW_ratio=[0.5 1 0];%set_PW_ratio=[1 0.5 0];

%% 
nHour=24;
nDays=365;
nTime=nDays*nHour;

%efficiency
%Battery
eff_total=0.8;
eff_c=eff_total^0.5;%利用可能な余剰電力を減らす
eff_d=1/eff_total^0.5;%埋めるべき残余需要を増やす　%→eff_c/eff_d=eff

%LDES
eff_total_LDES = 0.4;
eff_c_LDES = 0.8;
eff_d_LDES = 1/0.5;
%------------------
%% 制約条件設定
define_A_type_in = 1;%1:フリー、2：S1→S2のフロー＝０、3: ST2=0、4：ST1=0

%% コスト設定
%cg, c1e:tohoku2018_contour2の結果より2ケース設定（LCOE=10）
set_cg = [4.7];%yen/kWh singleLCOE
set_c1e = [0.5]*10^3;%yen/kWh

%cg, c1e, c1p, c2e, c2pin, c2pout
setCost= [NaN NaN 10000 10 10000 15000];
Is1p_yr = setCost(3);
Is2e_yr = setCost(4);
Is2pin_yr = setCost(5);
Is2pout_yr = setCost(6);

%Backupコスト設定%導入されないように、十分高価に設定
If = 10000000;
pf = 100000;

%% LP計算
rec_LCOE = NaN(n_reg, n_year, 3);%3ケース分(PV=0.5 1 0)
rec_x_set = NaN(n_reg*n_year, 7, 3);
rec_x_flow = cell(n_reg, n_year, 3);
rec_bottleneck = NaN(n_reg*n_year, 6, 3);

tic
for i_reg = set_i_reg
    for i_year = set_i_year
        for i_PVratio=1:length(set_PW_ratio)
        
            % プロファイル取得(Tohoku 2018ならregionNo=2,yearNo=3; Tokyo 2017ならregionNo=3,yearNo=2;）        
            regionNo = i_reg;
            yearNo = i_year;

            %規格化データ
            Demandunit=NormalizedData(:,1,regionNo,yearNo);
            PVunit=NormalizedData(:,2,regionNo,yearNo);
            Windunit=NormalizedData(:,3,regionNo,yearNo);
            sumD=sum(Demandunit);
            Gunit=(set_PW_ratio(i_PVratio)*PVunit + (1-set_PW_ratio(i_PVratio))*Windunit);%合成プロファイル（単位xgあたり）

            %コスト設定
            i_case = 1;
                LCOEg = set_cg(i_case)*sumD;
                Is1e_yr = set_c1e(i_case);
                Pin = [LCOEg, If, Is1e_yr, Is1p_yr, Is2e_yr, Is2pin_yr, Is2pout_yr, pf];            

                [x_set,x_flow, y, residue, S1toS2, sysLCOE, l_Dineq] = subfunc_B_LP_210723_xg_xs_xs2_system(Demandunit, Gunit, Pin, eff_c, eff_d, eff_c_LDES, eff_d_LDES, define_A_type_in);        

                rec_LCOE(i_reg, i_year, i_PVratio)=sysLCOE;
                rec_x_set((i_reg-1)*n_year + i_year,:, i_PVratio)=x_set;
                rec_x_flow{i_reg, i_year, i_PVratio} = x_flow;

                [out_x1e_x2e, term_ST2]=subfunc_bottleneck(x_set, Gunit, Demandunit, nTime, eff_c, eff_d, eff_c_LDES, eff_d_LDES);            
                rec_bottleneck((i_reg-1)*n_year + i_year,1:2+2*size(term_ST2,1), i_PVratio) = [out_x1e_x2e reshape(term_ST2',1,length(term_ST2(:)))];

                toc
                disp([i_reg i_year i_PVratio])            
        end
    end
end

%% result save
%save profile30_PV05_1_0.mat rec_LCOE rec_x_set rec_x_flow rec_bottleneck set_cg set_c1e setCost


%% ボトルネック期間特定
function [out_x1e_x2e, term_ST2]=subfunc_bottleneck(x_set, Gunit, Demandunit, nTime, eff_c, eff_d, eff_c_LDES, eff_d_LDES)
    %プロファイル集約（充放電効率による増幅なし）
    xg = x_set(1);
    r = xg*Gunit-Demandunit;
    [Rsum0, Gsum, ~, t_sum, ~] = M1_subfunction_residual_demand_sumup(r, Gunit, Demandunit);

    num_t_sum0 = t_sum(:,2)-t_sum(:,1)+1;%各集約期の期数
    num_t_sum0(num_t_sum0<0)=num_t_sum0(num_t_sum0<0)+nTime;%1に戻る集約期について期数の補正
    num_t_sum=zeros(length(num_t_sum0),2);
    num_t_sum(1:2:end-1,1)=num_t_sum0(1:2:end-1);%Plus期
    num_t_sum(2:2:end,2)=num_t_sum0(2:2:end);%Minus期
    clear num_t_sum0

    %RsumにST1の効率反映
    Rsum = eff_c*Rsum0.*(Rsum0>0) + eff_d*Rsum0.*(Rsum0<0);

    %CD matrix(充放電回数i(n)×初期集約プロファイルj(2n)）作成：
    [rec_xse, t_cd, rk_cell, CD, iMat]  = function_CDmatrix(Rsum, Gsum, t_sum, eff_d);    

    %LP最適解に対応するk（ダルマ落とし回数）
    k_pick = find(abs(rec_xse(:,1)-x_set(3)) == min(abs(rec_xse(:,1)-x_set(3))));
    %決定期間ST1
    tk = t_cd{k_pick};

    rk_pick0 = rk_cell{k_pick};
    rk_pick = rk_pick0.*(rk_pick0>0)*eff_c_LDES/eff_c + rk_pick0.*(rk_pick0<0)*eff_d_LDES/eff_d;
    nTsum=length(rk_pick);
    Hr=repmat(rk_pick,1,nTsum);%列ベクトルをn本並べてn×n正方行列化
    Gr=zeros(nTsum*2,nTsum);
    Gr(1:nTsum,:)=tril(Hr);
    Gr(nTsum+1:end,:)=triu(Hr,1);%G=[tril(H);triu(H,1)];%Hの下三角行列、上三角行列のつなぎ合わせ→各日を始点とする1年分のｒの生成
    Fr=cumsum(Gr);clear Hr Gr

    [minFr] = min(Fr(:));%xs2e解
    %x_set(5)+minFr%x_set(5)=minFrの確認

    i_min_all = find(abs(Fr(:)-minFr)<10^-8); 
    [iFr, jFr]=ind2sub(size(Fr),i_min_all);%Ｆrの最小値をとるインデックス取得（確認用）
    if iFr>size(Fr,2)
        iFr = iFr - size(Fr,2);
    end

    term_ST2 = [tk(jFr,1) tk(iFr,2)];%ST2のボトルネック期＝2つ存在する場合もあり（端点解のため）
    out_x1e_x2e = [rec_xse(k_pick,1) -minFr];%xgのダルマ落としで得られるx1e, x2e
end


%% 集約化のためのプロファイル前準備
function [Rin, Gin, Din] = profile_in_generator(xg, Gunit, Demandunit, eff_c, eff_d, nTime)
    r0 = xg * Gunit - Demandunit;
    r0(abs(r0)<10^-12) = 0;
    J0 = r0<0;%Q0に属するｊ（original) 
    %Q0 = sum( xg*Gunit(J0) )+sum( Demandunit(~J0) );%Q0：xse=0での再エネ利用量＝ J0期の発電量　＋　~J0期の需要量の合計

    eff_vector=ones(nTime,1)*eff_c;%rplusに対して減量の係数
    eff_vector(J0)=eff_d;%rminusに対して増量の係数                    
    Gin=eff_vector.*Gunit;%発電プロファイルに対して 同様の係数による増・減
    Din=eff_vector.*Demandunit;%需要プロファイルに対して 同様の係数による増・減
    clear eff_vector

    rplus=r0.* (~J0);
    rminus=r0.* J0;

    %充放電の入出力制約を考慮した余剰電力・残余需要プロファイル作成
    rplus(~J0) = eff_c * rplus(~J0);
    rminus(J0) = eff_d * rminus(J0);

    Rin = rplus + rminus;
end
%% 残余需要の集約化（＋－交互に登場. 第1期は＋）
function [Rsum, Gsum, Dsum, t_sum, Rdetail]=M1_subfunction_residual_demand_sumup(Rin, Gin, Din)

    %Rプラスorゼロ領域の各始点・終点の行番号
    t_Rplus_start=find( (Rin>=0 & [-1; Rin(1:end-1)]<0)==1 );
    t_Rplus_end=find( (Rin>=0 & [Rin(2:end); -1]<0)==1 );
    %Rマイナス領域の各始点・終点
    t_Rminus_start=find( (Rin<0 & [0; Rin(1:end-1)]>=0)==1 );
    t_Rminus_end=find( (Rin<0 & [Rin(2:end); 0]>=0)==1 );
            
    [Rsum0, Rdetail0]=function_profile_aggregation(Rin, t_Rplus_start, t_Rplus_end, t_Rminus_start, t_Rminus_end, Rin);%集約化ベクトル、各集約期ごとの元データ（並べ替え済み）
    Gsum0=function_profile_aggregation(Gin, t_Rplus_start, t_Rplus_end, t_Rminus_start, t_Rminus_end, Rin);
    Dsum0=function_profile_aggregation(Din, t_Rplus_start, t_Rplus_end, t_Rminus_start, t_Rminus_end, Rin);
    t_sum0 = function_t_aggregation(t_Rplus_start, t_Rplus_end, t_Rminus_start, t_Rminus_end, Rin, length(Rsum0));
    
    [j_start]=function_start_identification(Rsum0);%第1期とすべき期番号
    if j_start>1%残余需要の期の並べ替え
        Rsum = Rsum0([j_start:end 1:j_start-1]);
        Gsum = Gsum0([j_start:end 1:j_start-1]);
        Dsum = Dsum0([j_start:end 1:j_start-1]);
        t_sum = t_sum0([j_start:end 1:j_start-1], :);
        Rdetail =Rdetail0([j_start:end 1:j_start-1]);
    else
        Rsum = Rsum0;
        Gsum = Gsum0;
        Dsum = Dsum0;
        t_sum = t_sum0;
        Rdetail = Rdetail0;
    end
    
   

    %プロファイル集約に対応するオリジナルの期番号
    function [t_sum] = function_t_aggregation(t_Rplus_start, t_Rplus_end, t_Rminus_start, t_Rminus_end, Rin, nSum) 
        t_sum = NaN(nSum,2);
        if Rin(1)>=0
            if Rin(end)>=0%Vplusの最終期と1期の和→第1期
                t_sum(1:2:end-1,:) = [t_Rplus_start([end 2:end-1]) t_Rplus_end([1 2:end-1])];
                t_sum(2:2:end,:) = [t_Rminus_start t_Rminus_end];
            else                
                t_sum(1:2:end-1,:) = [t_Rplus_start t_Rplus_end];
                t_sum(2:2:end,:) = [t_Rminus_start t_Rminus_end];
            end
        else%Rin(1)<0
            if Rin(end)>=0%Vminusの第1期→最終期
                t_sum(1:2:end-1,:) = [t_Rplus_start t_Rplus_end];
                t_sum(2:2:end,:) = [t_Rminus_start([2:end 1]) t_Rminus_end([2:end 1])];
            else%Vminusの最終期と1期の和→最終期
                t_sum(1:2:end-1,:) = [t_Rplus_start t_Rplus_end]; 
                t_sum(2:2:end,:) = [t_Rminus_start([2:end-1 end]) t_Rminus_end([2:end-1 1])];
            end
        end            
        
    end

    %累積残余需要を最大にする期を第1期とする並べ替え
    function [j_start]=function_start_identification(Rsum)
        Rmat0 = repmat(Rsum, 1, length(Rsum));%Rsumをn列並べた正方行列
        Rmat = [tril(Rmat0,0); triu(Rmat0,1)];
        cumsumRmat_rplus_start = cumsum(Rmat(:, 1:2:end-1) ,1);%rplusが第1期となる列の累和
        [~, j_start0] = max(max(cumsumRmat_rplus_start,[],1));%最大値をとる列番号
        j_start = 2 * j_start0 -1;        
    end
end

%隣接するrマイナス・プラス領域の残余需要合計値の集計・記録
function [Vsum, Dcolumn] = function_profile_aggregation(V, t_Rplus_start, t_Rplus_end, t_Rminus_start, t_Rminus_end, Rin)        
    num_plus=length(t_Rplus_start);
    num_minus=length(t_Rminus_start);

    Vplus=zeros(num_plus,1);%余剰発電の集約
    Dplus=cell(num_plus,1);%余剰発電の集約期ごとのrtを降順並べ替えで記録
    for i=1:num_plus
        Vplus(i)=sum(V(t_Rplus_start(i):t_Rplus_end(i)));
        Dplus{i}=sort(V(t_Rplus_start(i):t_Rplus_end(i)), 'descend');
    end

    Vminus=zeros(num_minus,1);
    Dminus=cell(num_minus,1);%残余需要の集約期ごとのrtを降順並べ替えで記録
    for i=1:num_minus
       Vminus(i)=sum(V(t_Rminus_start(i):t_Rminus_end(i)));
       Dminus{i}=sort(V(t_Rminus_start(i):t_Rminus_end(i)), 'ascend');
    end           

    if Rin(1)>=0
        if Rin(end)>=0
            if length(Vminus)>1
                Vsum0=[[Vplus(end)+Vplus(1); Vplus(2:end-1)] Vminus];%Vplusの最終期と1期の和→第1期  
                D = [[{sort([Dplus{end}; Dplus{1}], 'descend')};Dplus(2:end-1)], Dminus]; 
            else
                Vsum0=[sum(Vplus) Vminus];%Vplusの最終期と1期の和→第1期  
                D = [{{sort([Dplus{end}; Dplus{1}], 'descend')}}, Dminus];                 
            end
        else
            Vsum0=[Vplus Vminus];                
            D = [Dplus, Dminus];
        end
    else%Rin(1)<0
        if Rin(end)>=0
            Vsum0=[Vplus Vminus([2:end 1])];%Vminusの第1期→最終期
            D = [Dplus, Dminus([2:end 1])];

        else
            Vsum0=[Vplus [Vminus(2:end-1);Vminus(end)+Vminus(1)]];%Vminusの最終期と1期の和→最終期                
            D = [Dplus, [Dminus(2:end-1); {sort([Dminus{end}; Dminus{1}], 'ascend')}]];
        end
    end            
    Vsum=reshape(Vsum0',[],1);%ベクトル化  
    Dcolumn = reshape(D', [], 1);
end

%% サブ関数
function [xg_rec_out, ab_rec_out]=sub_re100line_addtionalinfo(Demandunit,GENunit, rsum, xg, xg_plus)

    nTsum = length(Demandunit);
    
    %切片候補：Demandunit
    H0=repmat(Demandunit,1,nTsum);
    G0=zeros(nTsum*2,nTsum);
    G0(1:nTsum,:)=tril(H0);
    G0(nTsum+1:end,:)=triu(H0,1);%G=[tril(H);triu(H,1)];%Hの下三角行列、上三角行列のつなぎ合わせ→各日を始点とする1年分のｒの生成
    D=cumsum(G0);clear H0 G0 

    %傾き候補：GENunit
    H1=repmat(GENunit,1,nTsum);
    G1=zeros(nTsum*2,nTsum);
    G1(1:nTsum,:)=tril(H1);
    G1(nTsum+1:end,:)=triu(H1,1);%G=[tril(H);triu(H,1)];%Hの下三角行列、上三角行列のつなぎ合わせ→各日を始点とする1年分のｒの生成
    P=cumsum(G1);clear H1 G1
       
    
    % RE100におけるストレージ容量＝max(cumsum(rsum))    
    Hr=repmat(rsum,1,nTsum);%列ベクトルをn本並べてn×n正方行列化
    Gr=zeros(nTsum*2,nTsum);
    Gr(1:nTsum,:)=tril(Hr);
    Gr(nTsum+1:end,:)=triu(Hr,1);%G=[tril(H);triu(H,1)];%Hの下三角行列、上三角行列のつなぎ合わせ→各日を始点とする1年分のｒの生成
    Fr=cumsum(Gr);clear Hr Gr
    
    [minFr] = min(Fr(:));
    i_min_all = find(abs(Fr(:)-minFr)<10^-8);    
    
    %a,bの有効値抽出
    ab_rec=zeros(50*length(i_min_all),2);
    xgs_rec=zeros(50*length(i_min_all),2);   
    
    if length(i_min_all)>1
        disp(xg)%xgの探索では、i_min_allが複数は存在しない
    end
    
    for i = 1:length(i_min_all)
        i_min=i_min_all(i);
        [iFr, jFr]=ind2sub(size(Fr),i_min);%Ｆrの最小値をとるインデックス取得（確認用）
        
        %a,bの初期値抽出
        a0 = P(iFr,jFr);
        b0 = D(iFr,jFr);
        
        ab_rec(i,:)=[a0 b0];
        xgs_rec(i,:)=[xg -minFr];       
    end
    
    %a_vec, b_vec: a,bの選択候補 
    a_vec = P(:);
    b_vec = D(:);
    
    kk=length(i_min_all)+1;    
    xg_intersect = xg;
    while xg_intersect < xg_plus  && xg_intersect >= xg
        [max_ab]=max((a0-a_vec)./(b0-b_vec));
        I_max_all = find(abs((a0-a_vec)./(b0-b_vec) - max_ab)<10^-8);
        
        for ii=1:length(I_max_all)
            I_max=I_max_all(ii);
            b1=b_vec(I_max);
            a1=a_vec(I_max);        
            xg_intersect = (b0-b1)/(a0-a1);

            if xg_intersect < xg_plus && xg_intersect > xg  
                ab_rec(kk,:)=[a1 b1];
                xgs_rec(kk,:)=[xg_intersect -a1*xg_intersect+b1];
               
                kk=kk+1;
            else
                xg_intersect = 0;
            end            
        end
    end   
    ab_rec_out = ab_rec(1:kk-1,:);
    xg_rec_out = xgs_rec(1:kk-1,:);    
end

%% xse-Y関係 & CD matrix
function [rec_xse,t_cd, rk_cell, CD, iMat] = function_CDmatrix(Rsum, Gsum, t_sum, eff_d)
    nCalc = length(Rsum)/2;%集約プロファイルの山谷の数
        
    %record
    deltaXse=NaN(nCalc,1);%xseの増分
    deltaQdxg=NaN(nCalc,1);%dQdxg
    t_cd = cell(nCalc,1);%充放電の開始・終了期番号
    rk_cell = cell(nCalc,1);%r_k
    
    %初期プロファイル
    r_k_new=Rsum;    
    g_k_new=Gsum;    
    t_k_new=t_sum;
    
    %CD 充放電記録マトリクス
    CD=zeros(nCalc,length(Rsum));
    
    %Rsum_jが　ダルマ落とし後の集約プロファイルにおいて含まれる期番号の記録マトリクス
    iMat=zeros(nCalc,length(Rsum));
    
    
    %Mcalc0:正味需要対応マトリクスの更新前(n0とnk）
    Mcalc0=eye(nCalc*2);
           
    %端点となるxse計算＆集約プロファイルの更新
    for k=1:nCalc
        [dxse_k, dQdxg_k, t_kplus, Mcalc, t_min] = function_dalma_otoshi_matrix(r_k_new, g_k_new, t_k_new, Mcalc0);
        deltaXse(k)=dxse_k;
        deltaQdxg(k)=1/eff_d*dQdxg_k;
        t_cd{k}=t_k_new;
        rk_cell{k}=r_k_new;
        iMat(k,:)=[1:length(r_k_new)]*Mcalc0';
        
        q_k = Rsum' + sum(CD,1);%k-1での充放電後の正味需要(初期集約プロファイルベースのr_k_new）
        
        %Rsumへの配分法１:Rsumベースの残余比率
        temp_ratio = q_k ./ (Mcalc0*(q_k * Mcalc0)')';%Rsumベースの残余比率
        temp_ratio(isnan(temp_ratio))=0;%q_k=0で除する成分＝NaN→０に置換
        
        CD(k,1:2:end-1) = -dxse_k*temp_ratio(1:2:end-1);%充電
        CD(k,2:2:end) = dxse_k*temp_ratio(2:2:end);%放電        
                
        r_k_new=Mcalc'*(Rsum' + sum(CD,1))';        
        g_k_new=Mcalc'*Gsum;        
        t_k_new=t_kplus;
        
        Mcalc0 = Mcalc; 
    end
    
    %xse, dQdxse, Qの計算(出力変数)
    xse = cumsum(deltaXse);
    dQdxse=1/eff_d*[(nCalc:-1:1)';0];
    Q1 = cumsum( deltaXse.*dQdxse(1:end-1) );%残余需要への供給電力量    
    rec_xse = [[0;xse] [0;Q1]];
end

%ダルマ落としアルゴリズムによる、集約プロファイルに対する最小ストレージ量（kWh）の計算・プロファイルの更新
function [dxse_k, dQdxg_k, t_kplus, Mcalc, t_min] = function_dalma_otoshi_matrix(r_k, g_k, t_k, Mcalc0)
    %集約化残余需要プロファイルの絶対値最小値を探索＝xse_k(ST), 該当期番号 t_min
    [dxse_k, t_min]=min(abs(r_k));
        %※※複数成分がxse_kとなる場合も、min関数は最小のt_minのみ出力→
        %r_kplusにゼロの成分が発生する→dxse_k=0となる。

    %t_minでのg_kを場合分けして記録
    if r_k(t_min)>=0%min(r_k)は余剰発電
        dQdxg_k = g_k(t_min);
    else%残余需要がmin(r_k)
        dQdxg_k = - g_k(t_min);            
    end    
    
    %t_min-1:t_min+1の成分を統合化 →r_kplus, g_kplusに更新
    if t_min>1
        if t_min<length(r_k)% 1 < t_min < end        
            t_kplus = [t_k(1:t_min-2, :); [t_k(t_min-1, 1) t_k(t_min+1, 2)]; t_k(t_min+2:end,:)];            
            
            Mcalc = [Mcalc0(:,1:t_min-2) sum(Mcalc0(:,t_min-1:t_min+1),2) Mcalc0(:,t_min+2:end)];            
            
        else%t_min=end            
            Mcalc = [sum(Mcalc0(:, [end-1 end 1]),2) Mcalc0(:,2:end-2)];            
            
            if length(r_k)>2
                t_kplus = [[t_k(end-1, 1) t_k(1, 2)]; t_k(2:end-2,:)];                            
            else
                t_kplus = [t_k(1, 1) t_k(2, 2)];            
            end
        end
    else%t_min=1        
        t_kplus = [t_k(3:end-1,:); [t_k(end, 1) t_k(2, 2)]];       
        
        Mcalc = [Mcalc0(:,3:end-1) sum(Mcalc0(:, [end 1 2]),2)];                    
    end    
end

%% original B_LP_210705...
function [x_set,x_flow, y, residue, S1toS2, sysLCOE, l_Dineq]=subfunc_B_LP_210723_xg_xs_xs2_system(Demandunit, Gunit, Pin, eff_c1, eff_d1, eff_c2, eff_d2, define_A_type)
%価格条件------    
    Ig=Pin(1);
    If=Pin(2);
    Is1e=Pin(3);
    Is1p=Pin(4);
    Is2e=Pin(5);
    Is2pin=Pin(6);
    Is2pout=Pin(7);
    pf=Pin(8);


%% linprogオプション
%options = optimoptions('linprog','Algorithm','interior-point');
options = optimoptions('linprog','Algorithm','dual-simplex');
%options = optimoptions('linprog','Algorithm','dual-simplex', 'OptimalityTolerance',1*10^-8,'ConstraintTolerance',1*10^-6)
%% 最適化（充放電パターン）
%変数の並び　7 + 11*8760列
%設備(7)： xg xf xS1e xS1p xS2e xS2pin xS2pout 
%フロー：see notebook 21/07/05
%フロー[x1(G→(D)) x2(f→(D)) x3(G→(S1) S1に入る量はeff_c1掛け) x4(f→S1　S1に入る量はeff_c1掛け) x5(S1→(D) Dに入る量はeff_c1掛け) soc_1] * t 
%フロー2[x7(G→(S2) S2に入る量はeff_c2掛け) x8(f→S2　S2に入る量はeff_c2掛け) x9(S1-S2) eff_c1*eff_c2掛け  x10(S2→(D) Dに入る量はeff_c2掛け) soc_2] * t 

%制約条件行列（1期間）
%1）G
ag=[1 0 1 0 0 0 1 0 0 0 0];%(Gの発電量＞＝1*(G→D)　+　1*(G→B)）x1+x3+x7<=xgGt
%2）fire(backup)
af=[0 1 0 1 0 0 0 1 0 0 0];%x2 + x4 + x8 <= xf
%3）Demand
ad=[1 1 0 0 1/eff_d1 0 0 0 0 1/eff_d2 0];%(Demand　<=　1*(G→D)　+　1*(f→D) + eff_d*(B→D)）x1+x2+1/eff_d1*x5 + 1/eff_d2*x10 >= Dt

%4）S1制約(kWh)
b1=[0 0 eff_c1 eff_c1 -1 1 0 0 -1 0 0];%soc(t) eff_c*x3t + eff_c*x4t - x5t - x9t + soc(t) -soc(t+1) = 0
b2=[0 0 0 0 0 -1 0 0 0 0 0];%soc(t+1)
b3=[0 0 0 0 0 1 0 0 0 0 0];%soc(t)とxS1e or 0

%5）S2制約(kWh)
b4=[0 0 0 0 0 0 eff_c2 eff_c2 eff_c2/eff_d1 -1 1];%soc(t) eff_c*x3t + eff_c*x4t - x5t + eff_c2/eff_d1*x9t + soc(t) -soc(t+1) = 0
b5=[0 0 0 0 0 0 0 0 0 0 -1];%soc(t+1)
b6=[0 0 0 0 0 0 0 0 0 0 1];%soc(t)とxS2e or 0

%6) S1制約(kW)
s1_in=[0 0 1 1 0 0 0 0 0 0 0];%x3+x4 とxS1p or 0
s1_out=[0 0 0 0 1/eff_d1 0 0 0 1/eff_d1 0 0];%x5+x9 とxS1p or 0

%7) S2制約(kW)
s2_in=[0 0 0 0 0 0 1 1 1/eff_d1 0 0];%x7+x8+x9 とxS2pin or 0
s2_out=[0 0 0 0 0 0 0 0 0 1/eff_d2 0];%x10 とxS2pout or 0

%制約条件行列
n=length(Demandunit);

%speye(n)⇔eye(m)

    Ag=kron(speye(n),ag);
    Af=kron(speye(n),af);    

    B1=kron(speye(n),b1);
    B2_origin=kron(speye(n),b2);
    B2=[B2_origin(2:end,:); B2_origin(1,:)];%1期ずらし
    Beq_1=B1+B2;clear B1 B2 B2_origin

    B4=kron(speye(n),b4);
    B5_origin=kron(speye(n),b5);
    B5=[B5_origin(2:end,:); B5_origin(1,:)];%1期ずらし
    Beq_2=B4+B5;clear B4 B5 B5_origin
    
    D0=kron(speye(n),ad);

    %SOC制約　（Bu: -B+SOC<0, Bl:-SOC<0 ）
    Bu1=kron(speye(n),b3);
    Bl1=kron(speye(n),-b3);
    
    %SOC制約2　（Bu2: -B+SOC<0, Bl2:-SOC<0 ）
    Bu2=kron(speye(n),b6);
    Bl2=kron(speye(n),-b6);
    
    %第1ストレージkW制約
    S1_in = kron(speye(n), s1_in);
    S1_out = kron(speye(n), s1_out);
        
    %第2ストレージkW制約
    S2_in = kron(speye(n), s2_in);
    S2_out = kron(speye(n), s2_out);
    
    
    %制約行列統合--%特定x＝ゼロのオプション可能-----------
    %define_A_type
    switch define_A_type
        case 1
            %1）フリー
            A=sparse([-Gunit(1:n,1) sparse(n,6) Ag;...
                    sparse(n,1) -ones(n,1) sparse(n,5) Af]);%x9についての制約なし
        case 2        
            %2）x9＝0：S1→S2のフロー
            A=sparse([-Gunit(1:n,1) sparse(n,6) Ag;...
                    sparse(n,1) -ones(n,1) sparse(n,5) Af;
                    sparse(n,7) kron(speye(n),[zeros(1,8) 1 zeros(1,2)])]); %x9=0の制約
        case 3    
            %3）x10＝0→第2ストレージ＝０
            A=sparse([-Gunit(1:n,1) sparse(n,6) Ag;...
                    sparse(n,1) -ones(n,1) sparse(n,5) Af;
                    sparse(n,7) kron(speye(n),[zeros(1,9) 1 0])]); %x10=0の制約
        case 4
            %4）x5,x9＝0→第1ストレージ＝０
            A=sparse([-Gunit(1:n,1) sparse(n,6) Ag;...
                    sparse(n,1) -ones(n,1) sparse(n,5) Af;
                    sparse(n,7) kron(speye(n),[zeros(1,4) 1 0 0 0 1 0 0])]); %x5,x9=0の制約
                
        case 5
            %5）x3＝<B 第1ストレージの上限を設定
            A=sparse([-Gunit(1:n,1) sparse(n,6) Ag;...
                    sparse(n,1) -ones(n,1) sparse(n,5) Af;
                    sparse(1,2) 1 sparse(1,4)  zeros(1,11*n)]); %x3=<Bの制約
    end
   %--------------------------------------------------------
    D=sparse([sparse(n,7) -D0]); clear D0
   
    Bul1=sparse([sparse(n,2) -ones(n,1) sparse(n,4) Bu1;...
            sparse(n,7) Bl1]);clear Bu1 Bl1
    Bul2=sparse([sparse(n,4) -ones(n,1) sparse(n,2) Bu2;...
        sparse(n,7) Bl2]);clear Bu2 Bl2
    
    S_in_out = sparse([sparse(n,3) -ones(n,1) sparse(n,3) S1_in;...
                       sparse(n,3) -ones(n,1) sparse(n,3) S1_out;...
                       sparse(n,5) -ones(n,1) sparse(n,1) S2_in;...
                       sparse(n,6) -ones(n,1) S2_out]);clear S1_in S1_out S2_in S2_out               
    
    Beq=sparse([sparse(n,7) Beq_1;
                sparse(n,7) Beq_2]); clear Beq_1 Beq_2
    
    %目的関数係数
    f=sparse(zeros(1,7+n*11));    
    f(1:7)=[Ig If Is1e Is1p Is2e Is2pin Is2pout];
    f((7+2):11:((n-1)*11+7+2))=pf;
    f((7+4):11:((n-1)*11+7+4))=pf;
    f((7+8):11:((n-1)*11+7+8))=pf;
    
    %試行：
    %f((7+5):11:((n-1)*11+7+5))=0.01;
    %f((7+10):11:((n-1)*11+7+10))=0.011;
    %f((7+3):11:((n-1)*11+7+5))=0.01;
    %f((7+7):11:((n-1)*11+7+10))=0.011;
    %----------
    
    b=sparse([zeros(size(A,1),1);zeros(size(Bul1,1),1);zeros(size(Bul2,1),1);zeros(size(S_in_out,1),1); -Demandunit(1:n,1)]);
    beq=sparse([zeros(size(Beq,1),1)]);
    switch define_A_type
        case 5
            b(2*n+1)=15.55;
    end
    
   %[x, fval, exitflag, output, lambda] = linprog(f,[A;Bul1;Bul2;S_in_out;D],b,Beq,beq,sparse(zeros(n*11+7,1)),[],[],options);
   [x,fval,exitflag,output,lambda] = linprog_gurobi(f,[A;Bul1;Bul2;S_in_out;D],b,Beq,beq,sparse(zeros(n*11+7,1)),[],[]);
   
   x_set = x(1:7)';
   x_flow = [reshape(x(8:end),11,n)];
   
   %再エネ率 sum(x1+ eff1*x5 + eff2*x10)/sumD
   y = ( sum(x((7+1):11:((n-1)*11+7+1))) + 1/eff_d1*sum(x((7+5):11:((n-1)*11+7+5))) + 1/eff_d2*sum(x((7+10):11:((n-1)*11+7+10))))/sum(Demandunit);%再エネ率
   
   %余剰電力量 xg*sum(Gunit) - sum(x1 + x3 + x7)
   residue = x(1)*sum(Gunit) - sum(x((7+1):11:((n-1)*11+7+1))+x((7+3):11:((n-1)*11+7+3)) + x((7+7):11:((n-1)*11+7+7)));%余剰電力量
   
   S1toS2 = sum(x((7+9):11:((n-1)*11+7+9)));%ST1からST2へのフロー計
   
   sysLCOE = fval/sum(Demandunit);
   
   l_Dineq = lambda.ineqlin(end-n+1:end);
   
end
%% gurobiでLP
function [x,fval,exitflag,output,lambda] =linprog_gurobi(f,A,b,Aeq,beq,lb,ub,x0,ptions)
    % Build Gurobi model
    model.obj = full(f);
    model.A = [sparse(A); sparse(Aeq)]; % A must be sparse
    model.sense = [repmat('<',size(A,1),1); repmat('=',size(Aeq,1),1)];
    model.rhs = full([b(:); beq(:)]); % rhs must be dense
    if ~isempty(lb)
        model.lb = full(lb);
    else
        model.lb = -inf(size(model.A,2),1); % default lb for MATLAB is -inf
    end
    if ~isempty(ub)
        model.ub = full(ub);
    end

    % Extract relevant Gurobi parameters from (subset of) options
    params = struct();

%     if isfield(options,'Display') || isa(options,'optim.options.SolverOptions')
%         if any(strcmp(options.Display,{'off','none'}))
%             params.OutputFlag = 0;
%         end
%     end
% 
%     if isfield(options,'MaxTime') || isa(options,'optim.options.SolverOptions')
%         params.TimeLimit = options.MaxTime;
%     end

%計算過程表示OFF
params.OutputFlag = 0;

    % Solve model with Gurobi
    model.modelsense = 'min';
    result = gurobi(model,params);

    % Resolve model if status is INF_OR_UNBD
    if strcmp(result.status,'INF_OR_UNBD')
        params.DualReductions = 0;
        warning('Infeasible or unbounded, resolve without dual reductions to determine...');
        result = gurobi(model,params);
    end

    % Collect results
    x = [];
    output.message = result.status;
    output.constrviolation = [];

    if isfield(result,'x')
        x = result.x;
        if nargout > 3
            slack = model.A*x-model.rhs;
            violA = slack(1:size(A,1));
            violAeq = norm(slack((size(A,1)+1):end),inf);
            viollb = model.lb(:)-x;
            violub = 0;
            if isfield(model,'ub')
                violub = x-model.ub(:);
            end
            output.constrviolation = max([0; violA; violAeq; viollb; violub]);
        end
    end

    fval = [];

    if isfield(result,'objval')
        fval = result.objval;
    end

    if strcmp(result.status,'OPTIMAL')
        exitflag = 1; % converged to a solution
    elseif strcmp(result.status,'UNBOUNDED')
        exitflag = -3; % problem is unbounded
    elseif strcmp(result.status,'ITERATION_LIMIT')
        exitflag = 0; % maximum number of iterations reached
    else
        exitflag = -2; % no feasible point found
    end

    lambda.lower = [];
    lambda.upper = [];
    lambda.ineqlin = [];
    lambda.eqlin = [];

    if nargout > 4
        if isfield(result,'rc')
            lambda.lower = max(0,result.rc);
            lambda.upper = -min(0,result.rc);
        end
        if isfield(result,'pi')
            if ~isempty(A)
                lambda.ineqlin = -result.pi(1:size(A,1));
            end
            if ~isempty(Aeq)
                lambda.eqlin = -result.pi((size(A,1)+1):end);
            end
        end
    end

    if isempty(lambda.lower) && isempty(lambda.upper) && ...
            isempty(lambda.ineqlin) && isempty(lambda.eqlin)
        lambda = [];
    end

% Local Functions =========================================================
    
end
