how to make slack variables with only single index each but each comprises more than one element

I'm reproducing this paper and confused with the notation. There, z_{i} is a slack variable with single index. However, the definition of it: z_{i}: ith row vector in an m × N_{R} matrix of positive slack variables. My understanding is, each z_{i} has appropriate number of elements matching with (\zeta^{max})^T, i.e. if the dimension of (\zeta^{max})^T is 1 × 2 so z_{i} has to be 2 × 1 vector. How to write this slack variables in pyomo?

z_i, as you state, is a vector. So it must be double-indexed to represent the individual values. Anywhere you need to work with z_i as a vector in the equations, you'll need to manually break that down to representative expressions. For example the equations that have the dot products are just the sum products of z_i and another vector which can easily be written as a summation over the length of z_i.

Here's an example that mimics Eq 17 above. It is 1-based indexing like your example, which I find inconvenient, but matches the paper's structure.

Code

import pyomo.environ as pyo

M = 3  # rows
N = 4  # cols

model = pyo.ConcreteModel('vectors')

# SETS

model.M = pyo.Set(initialize=range(1, M+1))  # 1-based like example...
model.N = pyo.Set(initialize=range(1, N+1))

# VARS

model.z = pyo.Var(model.M, model.N, domain=pyo.Reals, doc='z vectors')


# from somewhere comes zeta...  
zeta = (1.0, 2.3, -5.8, -3.2)

# from somewhere comes d_i...
d = (1.0, 2.0, 3.0)


# CONSTRAINTS

# something similar to eqn 17 with a dot product...
# this is a "for each column i in M" type constraint

def c1(model, i):
    return - d[i-1] + sum(zeta[j-1] * model.z[i,j] for j in model.N) <= 0

model.eq17 = pyo.Constraint(model.M, rule=c1)

model.eq17.pprint()

Output:

eq17 : Size=3, Index=M, Active=True
    Key : Lower : Body                                                : Upper : Active
      1 :  -Inf : z[1,1] + 2.3*z[1,2] - 5.8*z[1,3] - 3.2*z[1,4] - 1.0 :   0.0 :   True
      2 :  -Inf : z[2,1] + 2.3*z[2,2] - 5.8*z[2,3] - 3.2*z[2,4] - 2.0 :   0.0 :   True
      3 :  -Inf : z[3,1] + 2.3*z[3,2] - 5.8*z[3,3] - 3.2*z[3,4] - 3.0 :   0.0 :   True
[Finished in 177ms]