Debugging cvxopt solve code

[Andrea Krijgsman]

Dear group,

See my code below for an optimization problem op  battery. When printing the number of charging/discharging time steps, I always get 5 and 5 regardless of how I set my parameters. It doesn't seem correct since there is no reason it is always the same number. Is there any mistake in my code that restricts that to max 10 switching times?

Thanks a lot!

Andrea

from cvxopt import matrix, solvers

from numpy import array, eye, hstack, ones, vstack, zeros

import numpy as np

from itertools import accumulate

import matplotlib.pyplot as plt

def generateG3(N):

    rb = np.zeros((N,N))

    lo = np.ones((N,N))*-1

    lb = np.zeros((N, N))

    for i in range(N):

        for j in range(i+1):

            lb[i,j] = -1

    ro = np.zeros((N, N))

    for i in range(N):

        for j in range(i+1):

            ro[i,j] = 1

   

    boven = np.hstack((lb,rb))

    onder = np.hstack((lo,ro))

    G3 = np.vstack((boven, onder))

    return G3

def generateG4(N):

    rb = np.zeros((N,N))

    lo = np.ones((N,N))

    lb = np.zeros((N, N))

    for i in range(N):

        for j in range(i+1):

            lb[i,j] = 1

    ro = np.zeros((N, N))

    for i in range(N):

        for j in range(i+1):

            ro[i,j] = -1

   

    boven = np.hstack((lb,rb))

    onder = np.hstack((lo,ro))

    G4 = np.vstack((boven, onder))

    return G4

def cvxopt_solve(N, solver=None):

    eff_c = 0.9

    eff_d = 1/ eff_c

    mean = 1.3

    var = 0.9

    p = np.random.normal(mean, np.sqrt(var), N)

    # charging decision variables

    c1 = [p[i]/eff_c for i in range(N)]

    # discharging decision variables

    c2 = [-p[i]*eff_d for i in range(N)]

    c = np.hstack((c1, c2))

    # cvxopt constraint format: G * x <= h

    # G1 constraint: x >= 0

    # there are n = 2N decision variables

    n = N*2

    G1 = -eye(n,n)

    h1  = zeros(n)

    # G2 constraint: x <= -0.1

    G2 = eye(n,n)

    x_min = 0.1

    h2  = ones(n)*x_min

    # G3 constraint: cumsum(x) >= 0 - SoC_start

    # -cumsum(x) <= -0 + SoC_start

    G3 = generateG3(N)

    min_cap = 0

    SoC_start = 0.5

    h3  = ones(n)*(-min_cap + SoC_start)

    # G4 constraint: cumsum(x) <= 1 - SoC_start

    G4 = generateG4(N)

    max_cap = 1

    h4  = ones(n)*(max_cap - SoC_start)

    # G5 constraint: sum(x) <= 0

    G5 = np.concatenate((np.ones(N), -np.ones(N)))

    h5  = zeros(1)

    # G6 constraint: sum(x) >= 0

    G6 = np.concatenate((-np.ones(N), np.ones(N)))

    h6  = zeros(1)

    # concantenate the constraints

    G = vstack([G1, G2, G3, G4, G5, G6])

    h = hstack([h1, h2, h3, h4, h5, h6])

    G = matrix(G)

    h = matrix(h)

    c = matrix(c)

    sol = solvers.lp(c, G, h, solver=solver)

    print('The price = ', [round(x, 2) for x in p])

    print('The optimal decision vector = ', [round(x, 2) for x in sol['x']])

    dSoC = [sol['x'][i] - sol['x'][i + n//2] for i in range(n//2)]

    print('The dSoC for each time step =', [round(x, 2) for x in dSoC])

    SoC = list(accumulate([SoC_start] + dSoC))

    print('The SoC per time step =', [round(x, 2) for x in SoC])

    CF = np.multiply(c , sol['x']).reshape(-1)

    print('Cashflow per decision step', [round(x, 2) for x in CF])

    CF2 = [CF[i] + CF[i + n//2] for i in range(n//2)]

    print('Cashflow per time step', [round(x, 2) for x in CF2])

    for i in range(N):

        print('p[{}] = {:.3f}, dSoC[{}] = {:.10f},  SoC[{}] = {:.3f}, CF[{}] = {:.3f}'.format(i, round(p[i], 3), i, round(dSoC[i],10),  i, round(SoC[i],3), i, round(CF2[i],6)))

    print('Check with objective = ', sum(CF2))

    print('The objective value = ', sol['primal objective'])

    print('The number of charging/discharging time steps = ', np.count_nonzero(dSoC))

    print('The number of charging time steps = ', sum(i > 0.01 for i in dSoC))

    print('The number of discharging time steps = ', sum(i < -0.01 for i in dSoC))

    fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)

    ax1.plot(SoC, color='blue')

    ax1.set_ylabel('SoC')

    ax1.set_title('State of Charge and Prices over Time')

    ax2.plot(p, color='green')

    ax2.set_xlabel('Time step')

    ax2.set_ylabel('Price')

    plt.show()    

cvxopt_solve(15, solver=None)