Create a matrix based on given constraints

I think you were on the right track with setting this up as an LP. It can be formulated as a MIP.

I haven't tinkered with any variety of inputs here, and I'm not sure you are guaranteed feasible results for all inputs with the constraints you have.

I penalized off-diagonal selection to encourage things on diagonal, and set up some "selection integrality" constraints to enforce block-selection.

Solves in about 1/10 of second...

# magic matrix

# Constraints
# 1) Column sum, CS, is: 300 <= CS <= 390
# 2) Row sum, RS, is equal to user-specified values, which are present in the E&T ticket column of the file
# 3) Number of non-zero values, NZV, in each column, should be: 0 < NZV <= 4
# 4) The NZV in the matrix should be: NZV >= 10
# 5) The pieces are stacked on top of each other. So, a the cell under a non-zero value cell is zero, than all cells underneath should have zeros.

import pyomo.environ as pyo

# user input
row_tots = [427.7, 12.2, 352.7, 58.3, 22.7, 31.9, 396.4, 29.4, 171.5, 474.5, 27.9, 200]
min_col_sum = 300
max_col_sum = 390
max_non_zero = 4
min_size = 10
bigM = max(row_tots)

m = pyo.ConcreteModel()

# SETS
m.I = pyo.Set(initialize=range(len(row_tots)))
m.I_not_first = pyo.Set(within=m.I, initialize=range(1, len(row_tots)))
m.J = pyo.Set(initialize=range(int(sum(row_tots)/min_col_sum)))

# PARAMS
m.row_tots = pyo.Param(m.I, initialize={k:v for k,v in enumerate(row_tots)})

# set up weights (penalties) based on distance from diagonal line
# between corners using indices as points and using distance-to-line formula
weights = { (i, j) : abs((len(m.I)-1)/(len(m.J)-1)*j - i) for i in m.I for j in m.J}
m.weight  = pyo.Param(m.I * m.J, initialize=weights)

# VARS
m.X = pyo.Var(m.I, m.J, domain=pyo.NonNegativeReals)
m.Y = pyo.Var(m.I, m.J, domain=pyo.Binary)          # selection indicator
m.UT = pyo.Var(m.I, m.J, domain=pyo.Binary)         # upper triangle of non-selects

# C1: col min sum
def col_sum_min(m, j):
    return sum(m.X[i, j] for i in m.I) >= min_col_sum
m.C1 = pyo.Constraint(m.J, rule=col_sum_min)

# C2: col max sum
def col_sum_max(m, j):
    return sum(m.X[i, j] for i in m.I) <= max_col_sum
m.C2 = pyo.Constraint(m.J, rule=col_sum_max)

# C3: row sum 
def row_sum(m, i):
    return sum(m.X[i, j] for j in m.J) == m.row_tots[i]
m.C3 = pyo.Constraint(m.I, rule=row_sum)

# C4: max nonzeros
def max_nz(m, j):
    return sum(m.Y[i, j] for i in m.I) <= max_non_zero
m.C4 = pyo.Constraint(m.J, rule=max_nz)


# selection variable enforcement
def selection_low(m, i, j):
    return min_size*m.Y[i, j] <= m.X[i, j]
m.C10 = pyo.Constraint(m.I, m.J, rule=selection_low)
def selection_high(m, i, j):
    return m.X[i, j] <= bigM*m.Y[i, j]
m.C11 = pyo.Constraint(m.I, m.J, rule=selection_high)

# continuously select blocks in columns.  Use markers for "upper triangle" to omit them

# a square may be selected if previous was, or if previous is in upper triangle
def continuous_selection(m, i, j):
    return m.Y[i, j] <= m.Y[i-1, j] + m.UT[i-1, j]
m.C13 = pyo.Constraint(m.I_not_first, m.J, rule=continuous_selection)
# enforce row-continuity in upper triangle
def upper_triangle_continuous_selection(m, i, j):
    return m.UT[i, j] <= m.UT[i-1, j]
m.C14 = pyo.Constraint(m.I_not_first, m.J, rule=upper_triangle_continuous_selection)
# enforce either-or for selection or membership in upper triangle
def either(m, i, j):
    return m.UT[i, j] + m.Y[i, j] <= 1
m.C15 = pyo.Constraint(m.I, m.J, rule=either)

# OBJ:  Minimze number of selected cells, penalize for off-diagonal selection
def objective(m):
    return sum(m.Y[i, j]*m.weight[i, j] for i in m.I for j in m.J)
#   return sum(sum(m.X[i,j] for j in m.J) - m.row_tots[i] for i in m.I) #+\
#           sum(m.Y[i,j]*m.weight[i,j] for i in m.I for j in m.J)
m.OBJ = pyo.Objective(rule=objective)
    

solver = pyo.SolverFactory('cbc')
results = solver.solve(m)

print(results)
for i in m.I:
    for j in m.J:
        print(f'{m.X[i,j].value : 3.1f}', end='\t')
    print()
print('\npenalty matrix check...')
for i in m.I:
    for j in m.J:
        print(f'{m.weight[i,j] : 3.1f}', end='\t')
    print()

Result

 300.0   127.7   0.0     0.0     0.0     0.0     0.0    
 0.0     12.2    0.0     0.0     0.0     0.0     0.0    
 0.0     165.6   187.1   0.0     0.0     0.0     0.0    
 0.0     0.0     58.3    0.0     0.0     0.0     0.0    
 0.0     0.0     22.7    0.0     0.0     0.0     0.0    
 0.0     0.0     31.9    0.0     0.0     0.0     0.0    
 0.0     0.0     0.0     300.0   96.4    0.0     0.0    
 0.0     0.0     0.0     0.0     29.4    0.0     0.0    
 0.0     0.0     0.0     0.0     171.5   0.0     0.0    
 0.0     0.0     0.0     0.0     10.0    390.0   74.5   
 0.0     0.0     0.0     0.0     0.0     0.0     27.9   
 0.0     0.0     0.0     0.0     0.0     0.0     200.0