[Algorithm] 도시화과정 시뮬레이션

 셀룰러 오토마타 (Cellular Automata) 는 유한한 격자 공간 내에서 세포(CA) 의 움직임을 구현한 알고리즘이다. 각 세포는 일정한 상태를 (예를 들면 생존/죽음 과 같은) 가지며, 시간의 흐름에 따라 일정한 규칙을 적용 받는다. 대표적 알고리즘 적용 사례로는 '생명게임' 이 있다.

 

 다음은 '생명게임' 을 구현한 코드와 결과이다.

import cellpylib as cpl


cellular_automaton = cpl.init_simple2d(60, 60)


cellular_automaton[:, [28,29,30,30], [30,31,29,31]] = 1 # Glider

cellular_automaton[:, [40,40,40], [15,16,17]] = 1 # Blinker

cellular_automaton[:, [18,18,19,20,21,21,21,21,20], [45,48,44,44,44,45,46,47,48]] = 1 # Light Weight Space Ship (LWSS)


# evolve the cellular automaton for 60 time steps
cellular_automaton = cpl.evolve2d(cellular_automaton, timesteps=60, neighbourhood='Moore',
                                  apply_rule=cpl.game_of_life_rule, memoize='recursive')

cpl.plot2d_animate(cellular_automaton)

생명게임 구현 결과

 여기서 행위자를 추가하여 각 행위자 간 상호작용을 관측하는 개념이 행위자 기반 모델 (Agent Based Model) 이다.

 

 

 개념이 중요한 것은 아니니 간략하게 설명만 하고 넘어간다. 중요한 것은

1. 행위자가 있다.
2. 행위자는 주어진 규칙에 의해 움직인다.

는 사실이다. 행위자 기반 모델에서 컨셉을 얻어서 도시화 과정을 시뮬레이션 하는 알고리즘을 만들어 본다.

 

 도시화 과정은 다음과 같다.

도시화 과정

 도시는 크게 [도시화 > 교외화 > 역도시화 > 재도시화] 과정을 거친다. 초기 도시가 형성될 때 인구가 집중되어 도시화가 일어나고, 높은 인구밀도로 인한 불경제로 인해 고소득층이 교외로 빠져나간다. 인구 유출이 점점 심화되며 중심도시 인구가 줄어들고 중심지가 쇠퇴하다가, 젠트리피케이션을 통해 중심지가 다시 부흥하고 유출된 인구가 다시 중심지로 들어온다.

 

 이제 이 일련의 과정에 참여하는 행위자와 이들의 행동을 설명할 수 있는 행동 규칙을 설정한다.

 

 우선 행위자는 소득에 따라 3 개 주체로 구분한다. [Low Income / Middle Income / High Income] 이들은 기본적으로 같은 행동 규칙을 갖지만, 세부적인 변수 설정에 차이를 두어 이동 양상을 구분했다. 이번 알고리즘에서는 각각 30 / 50 / 20 명으로 인구비율을 구성해서 랜덤으로 배치했다.

 

 행동 규칙은 크게 [Move / Stay ] 두 가지 선택지를 가진다. 인구밀도가 높아질수록 집적 이익과 불이익이 동시에 높아지지만 밀도가 일정 수준 이상이 되면 불이익이 더 커진다. 집적 이익과 불이익은 각각 무작위 선택의 가중치가 되어 의사 결정에 절대적 영향을 미친다. 인구밀도에 따른 집적 이익과 불이익을 나타낸 함수는 다음과 같다.

Advantage / Disadvantage Graph

파란 선은 집적 이익을, 빨간 선은 집적 불이익을 나타내는 그래프이다. 40을 기점으로 (총 인구수의 40%) 인구밀도가 너무 높다고 판단하여 집적 불이익이 더 높아지도록 함수를 설계했다. 여기에 소득에 따라 이동비용을 차등 적용했다. 저소득층은 이동이 자유롭지 않아 집적 불이익보다 집적 이익이 더 중요한 반면, 고소득층은 비교적 이동이 자유로우므로 집적 불이익에 더 민감하게 반응한다. 따라서 저소득층은 Advantage 함수에, 고소득층은 Disadvantage 함수에 각각 상수 10을 더해줬다.

 만약 Move 를 선택한 경우, 주변 8개 구역으로 이동한다. 이 때 각 구역의 인구밀도를 가중치로 두어 무작위 장소로 이동한다.

 

위 설정들을 적용한 알고리즘을 총 20번의 Time Step 동안 반복해 보았다.

초기 상태

Time Step 별 변화 양상 [ 0 > 6 > 14 > 18 ]

 초반엔 무작위로 흩어져 있는 주체들이 점차 한 곳으로 모여서 도시를 형성하는 것을 볼 수 있다. step 0 에서 가장 인구밀도가 높았던 곳이 점차 쇠퇴하고, 우상단에 위치한 도시가 점차 발달한다. 그 외 중소 도시들은 큰 변화 없이 제 자리를 지키는 것을 볼 수 있다. step 14 부터는 큰 변화 없이 기존 상태를 답보하는 양상을 확인했다.

 재밌는 점은, step6 부터 발달하기 시작한 도시의 인구 구성을 살펴보면, 고소득자의 비율이 높다. 반대로 쇠퇴하는 도시에 남아있는 인구구성은 저소득층과 중간소득층의 비율이 높다. 전형적인 쇠퇴도시의 교외화 현상이다.

 

 

 도시화 과정 같은 복잡한 과정을 비교적 간단한 수식으로만 나타내려고 하다 보니 대단한 결과를 볼 수는 없었다. 다만 시뮬레이션의 컨셉을 파악하고 추후 다른 프로젝트에 적용할 아이디어를 얻었다는 점, 그리고 여타 라이브러리 없이 numpy 로만 알고리즘을 구성했다는 점에서 의미를 찾을 수 있다.

 반대로 한계점은 성능이 그리 뛰어나지 않다는 점과 초기 무작위 분포에 따라 양상이 크게 달라진다는 점에서 부족함을 느꼈다. 이 점들을 개선하여 New Point 프로젝트에도 적용해 볼 생각이다.

 

#%%
import numpy as np
import random
import matplotlib.pyplot as plt
from tqdm import tqdm
import sympy as sy


#%%
# Create Field
field = np.zeros(shape=(1000, 1000))


# Functions
def create_agent(x, y, seed):
    random.seed
    agents = np.array([random.sample(range(1000), x), random.sample(range(1000), x), np.full((x), y)]).T
    return agents

def count_population(field, agents):
    field_count = field.copy()
    
    for a in agents:
        field_count[int(a[1]/100)*100 : int(a[1]/100)*100+100, int(a[0]/100)*100 : int(a[0]/100)*100+100] += 1
    return field_count

def population_func(pop):
    if pop < 5:
        advantage = 0
        disadvantage = 0
    elif pop <= 40:
        advantage = 2 * pop
        disadvantage = (1.1165)**pop - (1.1165)**5
    else:
        advantage = 3 * pop - 40
        disadvantage = (1.1165)**pop - (1.1165)**5
    
    return advantage, disadvantage

def neighbor_func(agent, neighbor):
    location = [int(agent[0]/100)*100, int(agent[1]/100)*100]
    
    if neighbor == 0:
        location[0] -= 100
        location[1] -= 100
    elif neighbor == 1:
        location[0] -= 100
    elif neighbor == 2:
        location[0] -= 100
        location[1] += 100
    elif neighbor == 3:
        location[1] -= 100
    elif neighbor == 4:
        location[1] += 100
    elif neighbor == 5:
        location[0] += 100
        location[1] -= 100
    elif neighbor == 6:
        location[0] += 100
    else:
        location[0] += 100
        location[1] += 100
        
    return location

def moving_func(field, agents):
    new_agents = []
    
    for a in agents:
        pop = field[a[0], a[1]]
        
        field_state = np.unique(field)
        
        if pop < 5:
            ms = 'move'
            
            dtn = random.choices(field_state, weights=field_state)
            dtn = np.where(field == dtn)
            dtn = np.concatenate([dtn[0].reshape(dtn[0].shape[0], 1), dtn[1].reshape(dtn[0].shape[0], 1)], axis=1)
            dtn = random.choices(dtn)[0]
            dtn = np.array([dtn[0], dtn[1], a[2]])
            
            new_agents.append(dtn)
            
        else:
            advantage, disadvantage = population_func(pop)
            
            if a[2] == 0:
                advantage += 10
            elif a[2] == 2:
                disadvantage += 10
            
            ms = random.choices(['move', 'stay'], weights=[advantage, disadvantage])
            
            if ms == 'move':
                a = agents[0]
                
                field_state = []
                field_state.append(field[a[0]-100, a[1]-100])
                field_state.append(field[a[0], a[1]-100])
                field_state.append(field[a[0]+100, a[1]-100])
                field_state.append(field[a[0]-100, a[1]])
                field_state.append(field[a[0]+100, a[1]])
                field_state.append(field[a[0]-100, a[1]+100])
                field_state.append(field[a[0], a[1]+100])
                field_state.append(field[a[0]+100, a[1]+100])
                                
                dtn = random.choices(range(0, 8), weights=field_state)[0]
                dtn = np.array(neighbor_func(a, dtn))
                x, y = random.choices(range(dtn[0], dtn[0]+100))[0], random.choices(range(dtn[1], dtn[1]+100))[0]
                dtn = np.array([x, y, a[2]])
                
                new_agents.append(dtn)
            else:
                new_agents.append(a)
      
    new_agents = np.array(new_agents)
    
    return new_agents


#%%
# Advantage / Disadvantage Function
x0 = np.linspace(5, 40, 200)
x1 = np.linspace(40, 50, 100)
x2 = np.linspace(5, 50, 300)
y0 = 2 * x0
y1 = 3 * x1 - 40
y2 = 1.1165 ** x2 - 1.1165 ** 5

plt.figure(figsize=(10, 10))
plt.grid(color='gray', alpha=.5)
plt.plot(x0, y0, c='Blue', label='y = 2x')
plt.plot(x1, y1, c='Blue', label='y = 3x - 40')
plt.plot(x2, y2, c='Red', label = 'y = (1.1165)^x - (1.1165)^5')
plt.legend(loc = 'upper right')
plt.savefig(r"C:\Users\PC\Desktop\주영준\진행중인 프로젝트\도시형성과정" + '/Advantage_Disadvantage Function.png', dpi=100)
plt.show()

del x0, x1, x2, y0, y1, y2


# Create Agents
agent_0 = create_agent(30, 0, 123)
agent_1 = create_agent(50, 1, 124)
agent_2 = create_agent(20, 2, 125)

agents = np.concatenate([agent_0, agent_1, agent_2])
field_count = count_population(field, agents)


# Initial State
plt.figure(figsize=(10, 10))
plt.title("Agent_Scatter_init")
plt.scatter(agent_0[:, 0], agent_0[:, 1], label='0')
plt.scatter(agent_1[:, 0], agent_1[:, 1], label='1')
plt.scatter(agent_2[:, 0], agent_2[:, 1], label='2', c="Red")
plt.legend(loc = 'upper right')
plt.imshow(field_count, cmap='Greens')
plt.savefig(r"C:\Users\PC\Desktop\주영준\진행중인 프로젝트\도시형성과정" + '/Agent_Scatter_init.png', dpi=100)
plt.show()


# Simulation
for t in tqdm(range(0, 20), leave=True):
    if t == 0:
        new_agents = moving_func(field_count, agents)
        new_field = count_population(field, new_agents)
    else:
        new_agents = moving_func(new_field, agents)
        new_field = count_population(field, new_agents)

    plt.figure(figsize=(10, 10))
    plt.title("Agent_Scatter_{0}".format(t))
    plt.scatter(new_agents[:30,0], new_agents[:30, 1], label='0')
    plt.scatter(new_agents[30:80,0], new_agents[30:80, 1], label='1')
    plt.scatter(new_agents[80:,0], new_agents[80:, 1], label='2', c="Red")
    plt.legend(loc = 'upper right')
    plt.imshow(new_field, cmap='Greens')
    plt.savefig(r"C:\Users\PC\Desktop\주영준\진행중인 프로젝트\도시형성과정" + '/Agent_Scatter_{0}.png'.format(t), dpi=100)
    plt.show()

del agent_0, agent_1, agent_2, t

소스코드 전체이다.