Code to assign persons by age to different households

I am generating a synthetic population for which number of households of each size and age composition of each household are known. I am trying to assign persons by age to each of these households.

The total number of persons in each age group (column sum) should sum to

Children   45196
Adult     148949
Senior     12195

while the total number of persons in households of each size (1-20) should sum to

1    2276
2    9366
3   23739
4   47636
5   42475
6   28338
7    3675
8    3728
9    3672
10   3830
11   3894
12   3792
13   3770
14   3710
15   3795
16   3648
17   3672
18   3744
19   3800
20   3780

I tried to code this in python as a set of linear equations. However, the presence of negative solutions does not help. The code is as below, how can it be modified to generate the total number of persons of each of the age groups in different households.

# Total Population
Population = 206340

# Number of Children, Adults and Seniors
Demography = np.array([ 45196, 148949,  12195])

# Number of households by size
Household_size_distribution = np.array([2276, 4683, 7913,11909,8495,4723,525,466,408,383,354,316,290,265,253,228,216,208,200,189])

# Probability that a person of a certain age group belongs to a household of a certain size
Age_Composition = np.array([[7.000e-04, 7.702e-01, 2.291e-01],
       [1.890e-02, 8.066e-01, 1.745e-01],
       [1.486e-01, 8.027e-01, 4.870e-02],
       [2.519e-01, 7.180e-01, 3.010e-02],
       [2.732e-01, 6.719e-01, 5.490e-02],
       [3.046e-01, 6.337e-01, 6.170e-02]])

# Store Age compositions
x = np.zeros((20,3),dtype=np.float)
x[:5,:] = Age_Composition[:5,:]
# Age composition same for households with more than 6 persons
x[5:20,:] = np.repeat(Age_Composition[5][np.newaxis,:], 15,0)

# Normalize the age compositions column-wise: Children, Adults and Seniors
y = np.zeros((20,3),dtype=np.float)
y[:,0] = x[:20,0]/np.sum(x[:20,0])
y[:,1] = x[:20,1]/np.sum(x[:20,1])
y[:,2] = x[:20,2]/np.sum(x[:20,2])

# Store Coefficients of 60 variables
w = np.zeros((23,60),dtype=np.float)
w[:20,:3] = x

z = np.zeros((20,3),dtype=np.float)
z[:,0] = y[:,0] 
w[20] = np.reshape(z,(60))

z = np.zeros((20,3),dtype=np.float)
z[:,1] = y[:,1] 
w[21] = np.reshape(z,(60))

z = np.zeros((20,3),dtype=np.float)
z[:,2] = y[:,2] 
w[22] = np.reshape(z,(60))

rollnumber = np.arange(0,60,step=3)
for i in range(20):
  w[i]=np.roll(w[i],rollnumber[i]) 

# Ax=B  
A = w
B = np.zeros((23,1),dtype=np.int)
B[:20,:] = new_Household_size_distribution[:,None]*np.arange(1,21)[:,None]
B[20:,:] = Demography[:,None]

# Solution to the set of linear equations
f=np.matmul(np.linalg.pinv(A),B)
f=np.reshape(f,60)

from Code to assign persons by age to different households