karmarkar.py

高等数学规划课程中需要用到的karmarkar标准型的转换以及投影尺度算法的实现，输出过程目标函数值和karmarkar函数值，迭代过程图表。##注意输入是一般形式标准型的参数集。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Date    : 2019-10-11 17:21:42
# @Author  : ${QU JINMING}
# @Email    : ${qjming97@163.com}


import os
import math
import sklearn
import numpy as np
import csv
import matplotlib.pyplot as plt


def trans_to_karmarkar（a b c m n）:

    # 约束矩阵A 等式右端b（列向量） 目标参数c（列向量） 约束数量m 变量数量n

    at = np.transpose（a）
    ct = np.transpose（c）
    bt = np.transpose（b）

    fu_b = -1 * b
    fu_c = -1 * c

    fu_at = -1*at
    fu_bt = -1*bt

    i_n = np.identity（n）
    zero_1 = np.zeros（（m 2*m+n））  # 第一行的0
    zero_2 = np.zeros（（n n））  # 第二行的0
    zero_3 = np.zeros（（1 n））  # 第三行的0
    # print（zero_3）
    zero_4 = np.array（[[0]]）

    e_total = np.ones（2*m+2*n）
    e_last = np.ones（（1 2*m+2*n+2））

    # print（b.shapec.shapezero_3.shape）

    b_gang = np.concatenate（（b c zero_4） axis=0）
    # print（b_gang.shape）
    fu_b_bang = -1*b_gang

    # 按行拼接再按列拼接
    A_plus1 = np.concatenate（（a zero_1） axis=1）
    # print（A_plus1.shape）
    A_plus2 = np.concatenate（（zero_2 at fu_at i_n） axis=1）
    # print（A_plus2.shape）
    A_plus3 = np.concatenate（（ct fu_bt bt zero_3） axis=1）
    # print（A_plus3.shape）
    A_plus = np.concatenate（（A_plus1 A_plus2 A_plus3） axis=0）
    # print（A_plus.shape）
    B_ae = b_gang - np.reshape（np.dot（A_plus e_total） （m+n+1 1））

    # print（B_ae.shape）

    B_1 = np.concatenate（（A_plus B_ae fu_b_bang） axis=1）
    # print（B_1.shape）
    B = np.concatenate（（B_1 e_last） axis=0）

    c_aim_1 = np.zeros（（2*m+2*n 1））
    temp = np.array（[[1] [0]]）
    c_aim = np.concatenate（（c_aim_1 temp） axis=0）

    return B c_aim


def karmarkar（B c_aim）:
    # 近似解的效果严重依赖于参数L 和 步长alpha
    # 收敛条件可采用L参数法 或比较目标函数值的变动比率
    # 调节两个参数的时候 alpha相当于learning rate 尽量小
    # L或者变动比率 L越小即绝对值越大收敛条件越苛刻 容易出现步长越界
    # 同理变动比率越小越容易越界 但越容易逼近目标最优值

    m = len（B）
    n = len（B[1]）
    A = B[:m-1]
    # print（“输入A“A）
    r = 1/math.sqrt（n*（n-1））
    # print（‘r‘r）
    # print（A）

    L = -22
    alpha = 1/32
    X_0 = np.ones（（n 1））/n
    # print（“x0“X_0）
    # X_0=np.array（[[0.6220][1/3][0.0447]]）
    Y_ORIGINAL = np.dot（np.transpose（c_aim） X_0）
    # print（“y00“Y_ORIGINAL）
    FLAG = Y_ORIGINAL*math.pow（2 L）  # 这个是教材给的收敛条件用的
    # print（“FLAG“FLAG）
    Y_each=[]
    X_each=[]
    print（“--------------迭代开始---------------\n“）
    while 1:
        X_1 = karmarkar_one_step（X_0 n A c_aim alpha r）
        Y_0 = np.dot（np.transpose（c_aim） X_0）
        Y_1 = np.dot（np.transpose（c_aim） X_1）

        assert Y_1 < Y_0  # 步长越界时该处报错

        X_0 = X_1
        Y_each.append（Y_1）
        X_each.append（X_1）
        print（“当前目标函数值：“ Y_1）

        bias = （Y_0-Y_1）/Y_1

        if bias < 0.0014:
            print（“--------------迭代结束---------------“）
            print（“最优近似变量值\n“ X_1 “\n最优近似 karmarkar目标值\n“ Y_1）
            break
        else:
            continue
    X_LAST = X_1