Решите систему N уравнений с N неизвестными, используя Джулию
У меня есть :
- набор из N локаций, которые могут быть местом работы или местом жительства
- вектор наблюдаемых рабочих L_i, с i в N
- вектор наблюдаемых жителей R_n, где n в N
- матрица расстояний между всеми парами мест жительства n и рабочих мест i
- эпсилон параметра формы
Установка N=3, эпсилон=5 и
d = [1 1.5 3 ; 1.5 1 1.5 ; 3 1.5 1] #distance matrix
L_i = [13 69 18] #vector of workers in each workplace
R_n = [27; 63; 10]
Я хочу найти вектор заработной платы (размер N), который решает эту систему N уравнений,
с л все рабочие места.
Нужно ли мне реализовывать итерационный алгоритм на векторах рабочих и заработной платы? Или можно напрямую решить эту систему?
Я пробовал это,
w_i = [1 ; 1 ; 1]
er=1
n =1
while er>1e-3
L_i = ( (w_i ./ d).^ϵ ) ./ sum( ( (w_i ./ d).^ϵ), dims=1) * R
er = maximum(abs.(L .- L_i))
w_i = 0.7.*w_i + 0.3.*w_i.*((L .- L_i) ./ L_i)
n = n+1
end
Ответы 1
Если L и R заданы (т. е. не зависят от w_i), вы должны настроить нелинейный поиск, чтобы получить (вектор) заработной платы из этого гравитационного уравнения (конечно, при условии нормализации одного w_i).
Вот минимальный пример. Я надеюсь, что это помогает.
# Call Packages
using Random, NLsolve, LinearAlgebra
# Set seeds
Random.seed!(1704)
# Variables and parameters
N = 10
R = rand(N)
L = rand(N) * 0.5
d = ones(N, N) .+ Symmetric(rand(N, N)) / 10.0
d[diagind(d)] .= 1.0
ε = -3.0
# Define objective function
function obj_fun(x, d, R, L, ε)
# Find shares
S_mat = (x ./ d).^ε
den = sum(S_mat, dims = 1)
s = S_mat ./ den
# Normalize last wage
x[end] = 1.0
# Define loss function
loss = L .- s * R
# Return
return loss
end
# Run optimization
x₀ = ones(N)
res = nlsolve(x -> obj_fun(x, d, R, L, ε), x₀, show_trace = true)
# Equilibrium vector of wages
w = res.zero
По закону Вальраса одно из N уравнений избыточно. Поэтому вы можете получить свой w_i только до нормализации. Я нормализовал последний, но в этом примере не имеет большого значения, какой вы выберете. См., например, обсуждение ниже уравнения (6) на странице 201 этой статьи: shorturl.at/yEFR5
@FabrizioLeone: il link che hai messo non funziona - ссылка на предоставленный вами документ не работает
Спасибо. Извините, попробуйте вот это
@FabrizioLeone: Grazie, но все равно Ошибка 404 ("Мы не можем найти страницу, которую вы ищете")...
Хорошо, еще одно испытание здесь. См. недатированную версию на веб-страницах авторов.
Какова цель x[end] = 1.0 ? Потому что, если я удалю его, я получу другой результат и извлеку результаты из своего итеративного алгоритма.