<commit_before>__author__ = 'vujanicr'
# when constructing the uncertainty vectors, it is important to know the order of the variables
# in the nominal problem in standard form; this is established in the __init__() method of the
# STN class, when filling the x_ijt, y_ijt and y_st dictionaries with cvx.Variable objects.

from STN import *

class robust_STN(object):
    """ robust_STN uses the nominal STN model, and provides
    - example methods to construct appropriate uncertainty sets,
    - a method to build the Robust Counterpart (RC)
    - a method to solve the RC
    - plotting
    """
    def __init__(self, stn=STN()):
        """
        :type stn: STN
        """
        self.stn = stn
        self.m = self.stn.m_eq + self.stn.m_ineq

        # Optimization Variables
        # ----------------------
        self.x = cvx.Bool(self.stn.n_x, 1)
        self.v = cvx.Variable(self.stn.n_y, 1)
        self.Y = cvx.Variable(self.stn.n_y, self.stn.n_x)  # recourse matrix
        # for the auxiliary vars Phi and Psi we work with dictionaries to have sparse representations
        self.Phi = {}
        self.Psi = {}
        # TODO BIG-M for Psi_bar
        self.Psi_bar = 10000*np.ones((self.m,1))

        # uncertainty sets are stored as matrices, within a list of length n_x
        # (number of bools in the problem).
        self.W = [0 for x in range(self.stn.n_x)]

    def x_ijt_index_to_std(self, x_ijt_ix):
        """ transforms the [i,j,t] index of the variable x_ijt in the index of the optimization variable
         x in the standard model.
        :param x_ijt_ix: [i,j,t] as a list containing the index of x_ijt
        :return: std_ix, the index of the x_ijt variable in vector of the model in standard form
        """
        i = x_ijt_ix[0]
        j = x_ijt_ix[1]
        t = x_ijt_ix[2]
        return i*(self.stn.J*self.stn.T) + j*self.stn.T + t

    def std_index_to_x_ijt(self, std_ix):
        """ transforms a given coordinate in the standard model into the appropriate [i,j,t] coordinate
        :param std_ix: the position within the variable x in the standard model
        :return: [i,j,t] as a list (ints)
        """
        i = int(std_ix/(self.stn.J*self.stn.T))
        remainder = std_ix%(self.stn.J*self.stn.T)
        j = int(remainder/(self.stn.T))
        t = int(remainder%(self.stn.T))
        return [i,j,t]

    def build_uncertainty_set_for_unit_delay(self, units=(0,), tasks=(0,), delay=1, from_t=0, to_t=None):
        # TODO fix doc, it's tuple; you need a comma at the end if it is one single unit (cannot iterate otherwise)
        """ Returns the uncertainty sets corresponding to possible delays of any task (in `tasks') executed
        on a unit (in `units'). The delays amount to up to `delay' time steps. Note: this is not equivalent to
        delays of exactly `delay' time steps, the ``up to'' is important. A version for `exact delays', which
        for instance can model the swapping of a certain operation from day to night is also possible.
        :param units: a (list of) unit(s) from STN.units (string); default='Heater'
        :param tasks: a (list of) task(s) from STN.tasks (string); default='Heat'
        :param delay: time delay, in steps (int); default=1
        :param from_t: delays may occur from time step ``from_t'' (int); default=0
        :param to_t: delays may occur only upt to time step ''to_t'' (int); default=STN.T
        :return: list of matrices [W_k], k=0,...,n_x-1
        """
        if to_t is None:
            to_t = self.stn.T-delay

        W = [0 for x in range(self.stn.n_x)]

        for i in units:
            for j in tasks:
                for t in range(from_t, to_t):
                    W_k = np.zeros((self.stn.n_x, delay+1))
                    # one column of the W_k matrix should always be 0, to encode that no event occurs;
                    # declaration left for clarity
                    W_k[:,0] = 0
                    # on to the other columns of W_k
                    for col in range(1, delay+1):
                        # index where we place the ``-1'' and ``+1''
                        minus_index = self.x_ijt_index_to_std([i,j,t])
                        plus_index = self.x_ijt_index_to_std([i,j,t+col])
                        # uncertainty vector w is -1 at (i,j,t), and +1 at (i,j,t+1)
                        w = np.zeros(self.stn.n_x)
                        w[minus_index] = -1
                        w[plus_index] = +1
                        W_k[:,col] = w

                    # assign matrix to the collection (list) of uncertainties
                    W[minus_index] = W_k
        # we return the value (and do not change self.W from within the method), since one may want
        # to combine several different W's before computing the RC
        return W

    def build_robust_counterpart(self):
        # fill in the required columns in Phi and Psi (can only be done after a W is constructed)
        k_ix = []  # list of indices with non-zero entries in W[k_ix]
        for k in range(self.stn.n_x):
            if np.any(self.W[k]):
                # TODO do I need to iterate or can I assign a vector?
                self.Phi[k] = cvx.Variable(self.m,1)
                self.Psi[k] = cvx.Variable(self.m,1)
                k_ix.append(k)

        # Aggregate Nominal Model
        # -----------------------
        A = np.vstack((self.stn.A_eq,self.stn.A_ineq))
        B = np.vstack((self.stn.B_eq,self.stn.B_ineq))
        # b = np.hstack((self.stn.b_eq.T,self.stn.b_ineq.T))
        # b = np.array([b])
        # b = b.T
        D = A

        # Constraints
        # -----------
        constraints = []
        for i, row in enumerate(self.stn.A_ineq):
             constraints.append( self.stn.A_ineq[i,:]*self.x + self.stn.B_ineq[i,:]*self.v +
                                 np.sum( [self.Phi[k][i] for k in k_ix] ) <= self.stn.b_ineq[i] )

        for i, row in enumerate(self.stn.A_eq):
             constraints.append( self.stn.A_eq[i,:]*self.x + self.stn.B_eq[i,:]*self.v +
                                 np.sum( [self.Phi[k][i+self.stn.m_ineq] for k in k_ix] ) == self.stn.b_eq[i] )

        for k in k_ix:
            constraints.append( np.ones((k_ix.__len__(),self.m))*np.diag(self.Psi[k]) >=
                                ((B*self.Y + D)*self.W[k]).T )
            constraints.append(self.Phi[k] >= 0)
            constraints.append(self.Phi[k] <= self.x[k]*self.Psi_bar)
            constraints.append(self.Psi[k] - self.Phi[k] >= 0)
            constraints.append(self.Psi[k] - self.Phi[k] <= (1-self.x[k])*self.Psi_bar)

        objective = cvx.Minimize(self.stn.c_x*self.x + self.stn.c_y*self.v)
        self.robust_model = cvx.Problem(objective, constraints)


if __name__ == '__main__':
    from STN import *
    stn = STN()
    stn.solve()
    rSTN = robust_STN(stn)
    rSTN.W = rSTN.build_uncertainty_set_for_unit_delay()
    rSTN.build_robust_counterpart()
    rSTN.robust_model.solve(verbose=True, solver='GUROBI')
<commit_msg>TODO: robust STN, encode non-anticipativity<commit_after>__author__ = 'vujanicr'
# when constructing the uncertainty vectors, it is important to know the order of the variables
# in the nominal problem in standard form; this is established in the __init__() method of the
# STN class, when filling the x_ijt, y_ijt and y_st dictionaries with cvx.Variable objects.

from STN import *
import scipy

class robust_STN(object):
    """ robust_STN uses the nominal STN model, and provides
    - example methods to construct appropriate uncertainty sets,
    - a method to build the Robust Counterpart (RC)
    - a method to solve the RC
    - plotting
    """
    def __init__(self, stn=STN()):
        """
        :type stn: STN
        """
        self.stn = stn
        self.m = self.stn.m_eq + self.stn.m_ineq

        # Optimization Variables
        # ----------------------
        self.x = cvx.Bool(self.stn.n_x, 1)
        self.v = cvx.Variable(self.stn.n_y, 1)
        self.Y = cvx.Variable(self.stn.n_y, self.stn.n_x)  # recourse matrix
        # for the auxiliary vars Phi and Psi we work with dictionaries to have sparse representations
        self.Phi = {}
        self.Psi = {}
        # TODO BIG-M for Psi_bar
        self.Psi_bar = 100000*np.ones((self.m,1))

        # uncertainty sets are stored as matrices, within a list of length n_x
        # (number of bools in the problem).
        self.W = [0 for x in range(self.stn.n_x)]

    def x_ijt_index_to_std(self, x_ijt_ix):
        """ transforms the [i,j,t] index of the variable x_ijt in the index of the optimization variable
         x in the standard model.
        :param x_ijt_ix: [i,j,t] as a list containing the index of x_ijt
        :return: std_ix, the index of the x_ijt variable in vector of the model in standard form
        """
        i = x_ijt_ix[0]
        j = x_ijt_ix[1]
        t = x_ijt_ix[2]
        return i*(self.stn.J*self.stn.T) + j*self.stn.T + t

    def std_index_to_x_ijt(self, std_ix):
        """ transforms a given coordinate in the standard model into the appropriate [i,j,t] coordinate
        :param std_ix: the position within the variable x in the standard model
        :return: [i,j,t] as a list (ints)
        """
        i = int(std_ix/(self.stn.J*self.stn.T))
        remainder = std_ix%(self.stn.J*self.stn.T)
        j = int(remainder/(self.stn.T))
        t = int(remainder%(self.stn.T))
        return [i,j,t]

    def build_uncertainty_set_for_unit_delay(self, units=(0,), tasks=(0,), delay=1, from_t=0, to_t=None):
        # TODO fix doc, it's tuple; you need a comma at the end if it is one single unit (cannot iterate otherwise)
        """ Returns the uncertainty sets corresponding to possible delays of any task (in `tasks') executed
        on a unit (in `units'). The delays amount to up to `delay' time steps. Note: this is not equivalent to
        delays of exactly `delay' time steps, the ``up to'' is important. A version for `exact delays', which
        for instance can model the swapping of a certain operation from day to night is also possible.
        :param units: a (list of) unit(s) from STN.units (string); default='Heater'
        :param tasks: a (list of) task(s) from STN.tasks (string); default='Heat'
        :param delay: time delay, in steps (int); default=1
        :param from_t: delays may occur from time step ``from_t'' (int); default=0
        :param to_t: delays may occur only upt to time step ''to_t'' (int); default=STN.T
        :return: list of matrices [W_k], k=0,...,n_x-1
        """
        if to_t is None:
            to_t = self.stn.T-delay-1

        W = [0 for x in range(self.stn.n_x)]

        for i in units:
            for j in tasks:
                for t in range(from_t, to_t):
                    W_k = np.zeros((self.stn.n_x, delay+1))
                    # one column of the W_k matrix should always be 0, to encode that no event occurs;
                    # declaration left for clarity
                    W_k[:,0] = 0
                    # on to the other columns of W_k
                    for col in range(1, delay+1):
                        # index where we place the ``-1'' and ``+1''
                        minus_index = self.x_ijt_index_to_std([i,j,t])
                        plus_index = self.x_ijt_index_to_std([i,j,t+col])
                        # uncertainty vector w is -1 at (i,j,t), and +1 at (i,j,t+1)
                        w = np.zeros(self.stn.n_x)
                        w[minus_index] = -1
                        w[plus_index] = +1
                        W_k[:,col] = w

                    # assign matrix to the collection (list) of uncertainties
                    W[minus_index] = W_k
        # we return the value (and do not change self.W from within the method), since one may want
        # to combine several different W's before computing the RC
        return W

    def build_robust_counterpart(self):
        # fill in the required columns in Phi and Psi (can only be done after a W is constructed)
        k_ix = []  # list of indices with non-zero entries in W[k_ix]
        for k in range(self.stn.n_x):
            if np.any(self.W[k]):
                self.Phi[k] = cvx.Variable(self.m,1)
                self.Psi[k] = cvx.Variable(self.m,1)
                k_ix.append(k)

        # Aggregate Nominal Model
        # -----------------------
        A = scipy.sparse.vstack((self.stn.A_eq,self.stn.A_ineq))
        B = scipy.sparse.vstack((self.stn.B_eq,self.stn.B_ineq))
        D = A

        # Constraints
        # -----------
        constraints = []
        for i, row in enumerate(self.stn.A_eq):
             constraints.append( self.stn.A_eq[i,:]*self.x + self.stn.B_eq[i,:]*self.v +
                                 np.sum( [self.Phi[k][i] for k in k_ix] ) == self.stn.b_eq[i] )
        for i, row in enumerate(self.stn.A_ineq):
             constraints.append( self.stn.A_ineq[i,:]*self.x + self.stn.B_ineq[i,:]*self.v +
                                 np.sum( [self.Phi[k][i+self.stn.m_eq] for k in k_ix] ) <= self.stn.b_ineq[i] )

        for k in k_ix:
            n_k = self.W[k].shape[1]
            constraints.append( np.ones((n_k,self.m))*cvx.diag(self.Psi[k]) >= ((B*self.Y + D)*self.W[k]).T )
            constraints.append(self.Phi[k] >= 0)
            constraints.append(self.Phi[k] <= self.x[k]*self.Psi_bar)
            constraints.append(self.Psi[k] - self.Phi[k] >= 0)
            constraints.append(self.Psi[k] - self.Phi[k] <= (1-self.x[k])*self.Psi_bar)
            # TODO: these following constraints help a lot computationally.
            #constraints.append( self.v + self.Y*self.W[k][:,1] <=  np.max(self.stn.V_max))
            #constraints.append( self.v + self.Y*self.W[k][:,1] >=  0)

        objective = cvx.Minimize(self.stn.c_x*self.x + self.stn.c_y*self.v)
        self.robust_model = cvx.Problem(objective, constraints)


if __name__ == '__main__':
    from STN import *
    stn = STN()
    stn.solve()
    rSTN = robust_STN(stn)
    rSTN.W = rSTN.build_uncertainty_set_for_unit_delay()
    rSTN.build_robust_counterpart()
    rSTN.robust_model.solve(verbose=True, solver='GUROBI')

{'fraction_non_alphanumeric': 0.12258700859518107, 'fraction_numerical': 0.005917993518388051, 'mean_word_length': 2.001691689574963, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1299470', 'n_tokens_mistral': 4874, 'n_tokens_neox': 4721, 'n_words': 1368}