Upload New File

Codigos Xpress/VarepsilonConstraintMaxMCLP.mos

+model VCMaxMCLP
+uses "mmxprs";
+uses "mmsystem"
+uses "mmive"
+
+(!
+	Declaraciones necesarias para la lectura del fichero
+!)
+declarations
+  m,n, p, r:integer
+  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_esp.txt" ! algo ms de 800 segundos la primera (es la ms importante)
+  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 Varepsilon-Constraint
+!)
+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
+	y: array(pdemanda) of mpvar
+	u: array(pdemanda) of mpvar
+	x: array(pservicio) of mpvar
+	psAbiertos : array(1..tam, 1..p) of integer
+	!p=3
+end-declarations
+
+
+(!
+	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 (MCLPR)
+!)
+!setparam("XPRS_MAXTIME",-100)
+obj := sum(j in pdemanda) u(j)
+!obj := sum(j in pdemanda) dem(j)*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(j in pservicio) x(j) is_binary
+forall(i in pdemanda) y(i) is_binary
+maximize(obj)
+
+varep := getobjval
+writeln("La solucion optima es\t", varep)
+
+
+n := 1
+
+! 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)
+
+res4 := sum(i in pdemanda) u(i) >= varep
+
+maximize(obj2) 
+
+Pn(n,2) := getobjval
+Pn(n,1) := varep
+
+writeln("La solucion del MCLP = ", Pn(n,1), "\t corresponde a maximizar MCLPR a ", Pn(n,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 (MCLPR)
+!)
+
+
+(!
+	Busquemos el comienzo de varepIter para no hacer calculo redundante: 
+		Esto lo hacemos sacando el maximo para MCLP y su correspondiente optimo de Pareto en MCLPR
+!)
+
+res4 := 0 ! nos aseguramos de no arrastrar una restriccion previa
+maximize(obj2)
+
+res4:= sum(i in pdemanda) y(i) >= getobjval 
+
+maximize(obj)
+
+varepIter := getobjval
+(! 
+	Fin Busquemos el comienzo de varepIter
+!)
+
+
+(!
+	Recorremos todo el espectro de entre varepIter y varep para obtener todas las soluciones de la frontera de pareto
+!)
+while(varepIter <= varep)do 
+	n := n+1
+	
+	! Encontrar la solucion optima para MCLP imponiendo como restriccion que el MCLPR sea varepIter
+	obj := sum(j in pdemanda) y(j)
+	!obj := sum(j in pdemanda) dem(j)*y(j)
+	
+	res4 := sum(i in pdemanda) u(i) >= varepIter
+	
+	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)
+	
+
+	! no queremos quedarnos con los puntos dominados 
+	if(  getobjval >= Pn(n-1,1)-0.1  ) then ! el menos una decima es necesario por culpa de que getobjval puede que sea un numero mal menos 10 a la menos el limite de XPRESS
+		n := n - 1
+		writeln("los datos del de anterior corresponden a un dominado")
+	end-if 
+	!writeln("")
+	!writeln("MCLP = ", getobjval, "\tMCLPR = ", varepIter, "-")
+	
+	! COMIENZO de una cosa nueva para quitar iteraciones de dominados (que no se hagan)
+	maxIterMCLP := getobjval
+	res4 := sum(i in pdemanda) y(i) >= maxIterMCLP
+	obj := sum(j in pdemanda) u(j)
+	
+	maximize(obj)	
+	
+	! FIN de una cosa nueva para quitar iteraciones de dominados (que no se hagan)
+	
+	
+	Pn(n,2) := getobjval
+	Pn(n,1) := maxIterMCLP
+	
+	varepIter := getobjval
+	writeln("MCLP = ", maxIterMCLP, "\tMCLPR = ", varepIter)
+	indR := 1
+	forall(indiceRellenar in pservicio | x(indiceRellenar).sol = 1) do 
+		psAbiertos(n, indR) := indiceRellenar
+		indR := indR +1 
+	end-do	
+	! incrementamos el valor de varepIter
+	varepIter := varepIter + 1
+
+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) writeln(Pn(ind2, 1),", ", Pn(ind2, 2)) 
+
+
+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), ", ")
+	end-do
+	writeln(" ")
+end-do
+
+
+
+
+!fichero de escritura de las soluciones: 
+!fopen("solucionVarepMCLP_large1.txt", F_OUTPUT)  ! Abre el archivo en modo escritura
+!fopen("solucionVarepMCLP_cyl.txt", F_OUTPUT)  ! Abre el archivo en modo escritura
+!fopen("solucionVarepMCLP_esp4.txt", F_OUTPUT)  ! Abre el archivo en modo escritura
+!fopen("solucionVarepMCLP_cyl_v2.txt", F_OUTPUT)  ! Abre el archivo en modo escritura
+fopen("solucionVarepMCLP_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
\ No newline at end of file