From 290be46343de2c4d09fc9173510bc61770d0584f Mon Sep 17 00:00:00 2001
From: tomcarr <tomas.carretero@estudiantes.uva.es>
Date: Sat, 29 Mar 2025 09:03:59 +0000
Subject: [PATCH] Upload New File
---
Codigos Xpress/NISE.mos | 472 ++++++++++++++++++++++++++++++++++++++++
1 file changed, 472 insertions(+)
create mode 100644 Codigos Xpress/NISE.mos
diff --git a/Codigos Xpress/NISE.mos b/Codigos Xpress/NISE.mos
new file mode 100644
index 0000000..5d44ae8
--- /dev/null
+++ b/Codigos Xpress/NISE.mos
@@ -0,0 +1,472 @@
+model NISE
+uses "mmxprs";
+uses "mmsystem"
+uses "mmive"
+
+(!
+ Declaraciones necesarias para la lectura del fichero
+!)
+declarations
+ m,n, p:integer
+ r: real
+ !archivo_datos = "instancesMOFLP/instances/large1.txt" ! formato_datos = 3
+ !archivo_datos = "municipios/archivo_salida.txt"
+ archivo_datos = "municipios/archivo_cyl.txt"
+ !archivo_datos = "municipios/archivo_esp2.txt"
+ formato_datos = 3
+
+
+end-declarations
+
+fopen(archivo_datos,F_INPUT)
+if(formato_datos = 1)then
+ read(n)
+ m:=n
+elif(formato_datos =2)then
+ read(m)
+ read(n)
+elif(formato_datos = 3)then
+ readln(m,n, p, r)
+elif(formato_datos = 4)then
+ readln(n)
+ m:=n
+end-if
+
+(!
+ Declaraciones necesarias para la lectura del fichero y el propio modelo NISE
+!)
+declarations
+ pdemanda= 1..n
+ pservicio = 1..m
+ cxS, cyS:array(pservicio)of real
+ cx, cy: array(pdemanda) of real
+ dist:array(pdemanda,pservicio)of real
+ !nc:array(pdemanda)of integer
+ dem:array(pdemanda)of real
+ a:array(pdemanda,pservicio)of integer
+ dc=r
+ tam = 1000
+ Pn : array(1..tam, 1..2) of real
+ Sn : array(1..tam, 1..2) of real
+ Snn : array(1..tam, 1..2) of real
+ Phis : array(1..tam) of real
+ Phisn : array(1..tam) of real
+ y: array(pdemanda) of mpvar
+ u: array(pdemanda) of mpvar
+ x: array(pservicio) of mpvar
+ psAbiertos : array(1..tam, 1..p) of integer
+ B : real
+ !p=3
+
+
+end-declarations
+
+(!
+ Funcion que dada la longitud de 3 lados de un triangulo, devuelve su area
+!)
+function area(a: real, b: real, c: real): real
+ declarations
+ s: real
+ resultado: real
+ end-declarations
+ s := (a + b + c) / 2
+ resultado := sqrt(s * (s - a) * (s - b) * (s - c))
+ returned := abs(resultado)
+end-function
+
+
+(!
+ Funcion que devuvelve la distancia entre dos puntos bidimensionales dados como sigue P1x, P1y, P2x, P2y
+!)
+function distEucl(x1: real, yp1: real, x2: real, y2: real): real
+ declarations
+ dist : real
+ end-declarations
+ dist := sqrt((x2 - x1)^2 + (y2 - yp1)^2)
+ returned := abs(dist)
+end-function
+
+(!
+ Terminamos de leer el fichero de datos y calculamos que puntos de servicio pueden dar servicio a que puntos demanda
+!)
+if(formato_datos = 1)then
+ forall(j in pservicio)read(j,cx(j),cy(j))
+ forall(i in pdemanda,j in pservicio)dist(i,j):=sqrt((cx(i)-cx(j))^2+(cy(i)-cy(j))^2)
+ forall(i in pdemanda,j in pservicio)do
+ if(dist(i,j)<=dc)then
+ a(i,j):=1
+ else
+ a(i,j):=0
+ end-if
+ end-do
+ forall(j in pservicio)costo(j):=1
+elif(formato_datos=2)then
+ forall(i in pdemanda)read(dem(i))
+ forall(i in pdemanda,j in pservicio)read(dist(i,j))
+ forall(i in pdemanda,j in pservicio)do
+ if(dist(i,j)<=dc)then
+ a(i,j):=1
+ else
+ a(i,j):=0
+ end-if
+ end-do
+ forall(j in pservicio)costo(j):=1
+elif(formato_datos = 3)then
+ forall(j in pservicio)readln(cxS(j),cyS(j))
+
+ forall(i in pdemanda)readln(cx(i),cy(i),dem(i))
+
+ forall(i in pdemanda,j in pservicio)dist(i,j):=sqrt((cx(i)-cxS(j))^2+(cy(i)-cyS(j))^2)
+
+ forall(i in pdemanda,j in pservicio)do
+ if(dist(i,j)<=r)then a(i,j):=1
+ else a(i,j):=0
+ end-if
+ end-do
+
+ elif(formato_datos = 4)then
+ forall(i in pdemanda)readln(cx(i),cy(i),dem(i))
+ forall(i in pdemanda,j in pservicio)do
+ dist(i,j):=round(sqrt((cx(i)-cx(j))^2+(cy(i)-cy(j))^2))
+ if(dist(i,j)<=dc)then a(i,j):=1
+ else a(i,j):=0
+ end-if
+ end-do
+end-if
+fclose(F_INPUT)
+
+
+
+(!
+ Comienzo maximizando el primer criterio (MCLP)
+!)
+!setparam("XPRS_MAXTIME",-100)
+obj := sum(j in pdemanda) y(j)
+!obj := sum(j in pdemanda) dem(j)*y(j)
+
+forall(i in pdemanda)
+ res1(i):= sum(j in pservicio)x(j)*a(i,j)>= y(i)
+
+forall(i in pdemanda)
+ res2(i):= sum(j in pservicio)x(j)*a(i,j) >= 2*u(i)
+
+res3 := sum(j in pservicio) x(j) = p
+
+forall(i in pdemanda) y(i) is_binary
+forall(j in pservicio) x(j) is_binary
+forall(i in pdemanda) u(i) is_binary
+
+maximize(obj)
+Pn(1,1) := getobjval
+
+! Busco para el otro criterio (MCLPR) la solucion optima, si pongo como restriccion que el MCLP sea el maximo anteiror
+obj2 := sum(i in pdemanda) u(i)
+
+resEliminar := sum(i in pdemanda) y(i) >= Pn(1,1)
+
+maximize(obj2)
+Pn(1,2) := getobjval
+writeln("")
+writeln("Max MCLP = ", Pn(1,1), "\t con MCLPR = ", Pn(1,2))
+! guardo los puntos abiertos para esta solcuion de la frontera de pareto
+indR := 1
+forall(indiceRellenar in pservicio | x(indiceRellenar).sol = 1) do
+ psAbiertos(1, indR) := indiceRellenar
+ indR := indR +1
+end-do
+(!
+ Fin maximizar el primer criterio (MCLP)
+!)
+
+
+
+(!
+ Paso a maximizar el segudno criterio (MCLPR)
+!)
+
+obj := sum(j in pdemanda) u(j)
+!obj := sum(j in pdemanda) dem(j)*u(j)
+
+resEliminar := sum(i in pdemanda) u(i) >= 0
+
+maximize(obj)
+
+Pn(2,2) := getobjval
+! Busco para el otro criterio (MCLP) la solucion optima, si pongo como restriccion que el MCLPR sea el maximo anteiror
+obj2 := sum(i in pdemanda) y(i)
+
+resEliminar := sum(i in pdemanda) u(i) >= Pn(2,2)
+
+maximize(obj2)
+Pn(2,1) := getobjval
+writeln("")
+writeln("MCLP = ", Pn(2,1), "\t correspondiente a Max MCLPR = ", Pn(2,2))
+indR := 1
+forall(indiceRellenar in pservicio | x(indiceRellenar).sol = 1) do
+ psAbiertos(2, indR) := indiceRellenar
+ indR := indR +1
+end-do
+(!
+ Fin de maximizar el segundo criterio
+ !)
+
+
+
+! Limpio las restricciones de tipo Eliminar
+resEliminar := 0
+
+
+
+(!
+ Guardo las dos primeras solcuiones de pareto previamente obtenidas de manera ordenada en S,
+ asi como calculo el punto intermedio del criterio NISE y pongo n = 2 (2 soluciones)
+!)
+imax := 2
+Sn(1,1) := Pn(2,1)
+Sn(1,2) := Pn(2,2)
+Sn(2,1) := Pn(1,1)
+Sn(2,2) := Pn(1,2)
+interX := Pn(1,1)
+interY := Pn(2,2)
+n := 2
+
+! Calculo la distiancia entre las dos soluciones que tenemos
+distP1P2 := distEucl( Pn(1,1),
+ Pn(1,2),
+ Pn(2,1),
+ Pn(2,2) )
+
+! Calculo el valor del criterio NISE para las dos soluciones
+Phis(1) := (area( (interX-Pn(2,1)),
+ (interY-Pn(1,2)),
+ distP1P2 )*2)/distP1P2
+
+Phi := Phis(1)
+
+! Establezco el error maximo permitido, criterio de finalizacion del bucle.
+T := Phi *0.008
+!T := Phi *0.01
+writeln(T)
+
+
+(!
+ Se itera hasta que todas las soluciones ordenadas consecutivas tienen un criterio menor que T,
+ por lo que en cada iteracion se podra annadir una nueva solucion de la frontera de pareto de Pareto.
+!)
+while(Phi <> 0) do
+ ! Obtego el indice el Phi maximo porque de esas dos soluciones sale la proxima intermedia.
+ maximo := 0.0
+ forall(iden in 1..tam)do
+ if(Phis(iden) >= maximo)then
+ imax := iden
+ maximo := Phis(iden)
+ end-if
+ end-do
+ ! Solucion optima NISE para el imax anterior
+ obj := abs(Sn(imax,2) - Sn((imax+1),2))*(sum(j in pdemanda)y(j)) + abs(Sn((imax+1),1)-Sn(imax,1))*(sum(j in pdemanda)u(j))
+
+ forall(i in pdemanda)
+ res1(i):= sum(j in pservicio)x(j)*a(i,j) >= 2*u(i)
+
+ forall(i in pdemanda)
+ res2(i):= sum(j in pservicio)x(j)*a(i,j) >= y(i)
+
+ res3 := sum(j in pservicio) x(j) = p
+
+ forall(i in pdemanda) u(i) is_binary
+ forall(i in pdemanda) y(i) is_binary
+ forall(j in pservicio) x(j) is_binary
+
+ maximize(obj)
+
+ ! Obtengo el valor de la solucion optima para MCLP y MCLPR
+ Syi := 0
+ forall(i in pdemanda) do
+ iter :=0
+ forall(j in pservicio | a(i,j) = 1 and x(j).sol > 0.999)do
+ iter := iter + 1
+ end-do
+ if(iter >= 2)then
+ Syi := Syi + 1
+ end-if
+ end-do
+
+ Sxi := 0
+ forall(ii in pdemanda) do
+ iter :=0
+ forall(j in pservicio | a(ii,j) = 1 and x(j).sol > 0.999)do
+ iter := iter + 1
+ end-do
+ if(iter >= 1)then
+ Sxi := Sxi + 1
+ end-if
+ end-do
+
+
+ ! Calculo B, si mayor que la sol optima anteiror entonces entre estos puntos que nos da imax, no se pueden encontrar mas soluciones del conjutno de pareto, de otro modo se annade la solucion de esta iteracion
+ B := abs(Sn(imax,2) - Sn((imax+1),2))*Sn(imax,1)
+ B := B + abs(Sn((imax+1),1)-Sn(imax,1))*Sn(imax,2) + 0.1 ! el mas una decima, hay que ponerlo porque si es igual se mete aunque no debiera, eso es culpa de XPRESS no mia, en algunos PC va bien otros no.
+ !writeln("\t", Sn((imax),1),"\t", Sn((imax),2),"\t", Sn((imax+1),1),"\t", Sn((imax+1),2))
+ writeln(getobjval, "\t", Sxi, "\t", Syi, "\t B = ",B)
+
+ if( ( Sn((imax),1) = Sn((imax+1),1) and Sn((imax),2) = Sn((imax+1),2) ) or ( Sn((imax),1) = Sxi and Sn((imax),2) = Syi ) or ( Sxi = Sn((imax+1),1) and Syi = Sn((imax+1),2) ) )then
+ !writeln("Es por esto")
+ Phis(imax) := 0
+ else if( getobjval > B ) then
+ ! incrementamos el numero de soluciones en la frontera de pareto de pareto
+ n := n + 1
+
+ Pn(n,1) := Sxi
+ Pn(n, 2) := Syi
+
+
+ ! para la solucion nueva guardamos que puntos de servicio abre
+ indR := 1
+ forall(indiceRellenar in pservicio | x(indiceRellenar).sol = 1) do
+ psAbiertos(n, indR) := indiceRellenar
+ indR := indR +1
+ end-do
+
+ ! Actualizamos segun corresponde los vectores Phi y matriz S
+ iterador := 1
+ while(iterador <= n)do
+ if(iterador <= imax)then
+ Snn(iterador, 1) := Sn(iterador, 1)
+ Snn(iterador, 2) := Sn(iterador, 2)
+ elif(iterador = (imax+1)) then
+ Snn(iterador, 1) := Sxi
+ Snn(iterador, 2) := Syi
+ else
+ Snn(iterador, 1) := Sn((iterador-1), 1)
+ Snn(iterador, 2) := Sn((iterador-1), 2)
+ end-if
+ iterador := iterador + 1
+ end-do
+
+ Sn := Snn
+
+ iterador := 1
+ while(iterador < (n))do
+ if(iterador < imax)then
+ Phisn(iterador) := Phis(iterador)
+
+ elif(iterador = imax) then
+ interX := Sn((imax+1),1)
+ interY := Sn(imax,2)
+
+ distSiSimas1 := distEucl( Sn((imax),1),
+ Sn((imax),2),
+ Sn((imax+1),1),
+ Sn((imax+1),2))
+
+ Phisn(iterador) := area((interX-Sn((imax),1)),
+ (interY-Sn((imax+1),2)),
+ distSiSimas1)*2/distSiSimas1
+
+ elif(iterador = (imax+1)) then
+ interX := Sn((imax+2),1)
+ interY := Sn((imax+1),2)
+
+ distSiSimas2 := distEucl( Sn((imax+1),1),
+ Sn((imax+1),2),
+ Sn((imax+2),1),
+ Sn((imax+2),2))
+
+ Phisn(iterador) := area((interX-Sn((imax+1),1)),
+ (interY-Sn((imax+2),2)),
+ distSiSimas2 )*2/distSiSimas2
+
+ else
+ Phisn(iterador) := Phis(iterador-1)
+ end-if
+ iterador := iterador + 1
+ end-do
+
+ Phis := Phisn
+
+
+ else
+ Phis(imax) := 0
+ end-if
+ end-if
+
+ ! Se presupone que se termina en todas las iteraciones, pero con que
+ ! un phi sea mayor que T, se sigue buscando soluciones
+ Phi := 0
+ forall(idenTerminar in 1..(n-1)) do
+ if(Phis(idenTerminar) > T)then
+ Phi := 1
+ end-if
+ end-do
+ !writeln("la lista de phis queda: \t", Phis)
+ !forall( ind2 in 1..n)do
+ ! writeln(Sn(ind2, 1),", ", Sn(ind2, 2))
+ !end-do
+end-do
+
+
+! Dibujar los pintos de servicio y las asignaciones
+id1:=IVEaddplot("puntos", IVE_RED)
+id2:=IVEaddplot("lineas", IVE_BLACK)
+id3:= IVEaddplot("mejores", IVE_BLUE)
+!forall(i in pdemanda) IVEdrawpoint(id1, cx(i), cy(i))! Dibujamos los puntos
+!forall(i in pdemanda, j in pservicio | y(i).sol > 0.999 and x(j).sol > 0.999 and a(i,j) = 1 ) IVEdrawline(id2, cx(i), cy(i), cx(j), cy(j)) ! Dibujamos las lineas
+
+
+forall(i in 1..n)do
+ IVEdrawpoint(id3, Pn(i,1), Pn(i,2))
+end-do
+!IVEdrawpoint(id3, Pn(1,1), Pn(1,2))
+!IVEdrawpoint(id3, Pn(2,1), Pn(2,2))
+IVEdrawpoint(id1, 0, 0)
+
+writeln("\n\n")
+
+writeln("El conjunto de soluciones factibles es \n")
+
+forall( ind2 in 1..n)do
+ writeln(Pn(ind2, 1),", ", Pn(ind2, 2))
+ !fprintln(f, Pn(ind2, 1),", ", Pn(ind2, 2))
+end-do
+
+writeln("\n\nLos puntoes de servicio que se habren son: \n")
+
+
+forall( indImprimir in 1..n)do
+ forall( ind2 in 1..p) do
+ write( psAbiertos(indImprimir,ind2), ", ")
+ !fprint(f, psAbiertos(indImprimir,ind2), ", ")
+ end-do
+ writeln(" ")
+ !fprintln(" ")
+end-do
+
+
+!fichero de escritura de las soluciones:
+!fopen("solucionNISE_large1.txt", F_OUTPUT) ! Abre el archivo en modo escritura
+!fopen("solucionNISE_cyl.txt", F_OUTPUT) ! Abre el archivo en modo escritura
+!fopen("solucionNISE_esp.txt", F_OUTPUT) ! Abre el archivo en modo escritura
+fopen("solucionNISE_esp_Borrar.txt", F_OUTPUT) ! Abre el archivo en modo escritura
+writeln(cx)
+writeln(cy)
+writeln(cxS)
+writeln(cyS)
+writeln(n)
+writeln(r)
+writeln(p)
+
+forall( ind2 in 1..n)do
+ writeln(Pn(ind2, 1),", ", Pn(ind2, 2))
+end-do
+
+
+forall( indImprimir in 1..n)do
+ forall( ind2 in 1..p) do
+ write( psAbiertos(indImprimir,ind2), " ")
+ end-do
+ writeln(" ")
+end-do
+
+fclose(F_OUTPUT) ! Cierra el archivo para guardar los cambios
+
+end-model
--
GitLab